Rの行列と配列の扱い方のざっとした説明
CやFortranでは行列といっても多次元配列があるだけで、行列としての操作はサポートされていないのですが、Rでは行列の演算が標準で提供されています。
1 行列の作成
行列の成分となるベクトルと、行数と列数のどちらかがあれば行列をつくれます。両方指定しても問題ないです。
matrix(1:10, 2) [,1] [,2] [,3] [,4] [,5]
[1,] 1 3 5 7 9
[2,] 2 4 6 8 10matrix(1:10, ncol = 5) [,1] [,2] [,3] [,4] [,5]
[1,] 1 3 5 7 9
[2,] 2 4 6 8 10ベクトルの要素は列ごとに割り当てていきますが、byrow = TRUEをつけると行ごとに使っていきます。
matrix(1:10, 2, byrow = TRUE) [,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10作成時に名前をつけることもできます。行名だけ、列名だけをつけることもできます。
(m <- matrix(1:4, 2, dimnames = list(c("a", "b"), c("A", "B")))) A B
a 1 3
b 2 4matrix(1:4, 2, dimnames = list(c("a", "b"), NULL)) [,1] [,2]
a 1 3
b 2 4matrix(1:4, 2, dimnames = list(NULL, c("A", "B"))) A B
[1,] 1 3
[2,] 2 41.1 対角行列
対角行列の作成はdiagと言う関数を使うと便利です。
diag(3) [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1引数に数字ひとつを入れると、その大きさの単位行列になりますが、複数にすると対角成分として使います。
diag(1:3) [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 31.2 行と列の追加
行の追加はrbind、列の追加はcbindで行います。
rbind(m, -1:-2) A B
a 1 3
b 2 4
-1 -2cbind(m, -1:-2) A B
a 1 3 -1
b 2 4 -22 行列の情報
行列は、要素数、行数、列数、行名、列名の情報を取ることができます。
length(m)[1] 4nrow(m)[1] 2ncol(m)[1] 2rownames(m)[1] "a" "b"colnames(m)[1] "A" "B"attributes(m)$dim
[1] 2 2
$dimnames
$dimnames[[1]]
[1] "a" "b"
$dimnames[[2]]
[1] "A" "B"2.1 行列であるかの判別
is.*が用意されています。
is.matrix(m)[1] TRUEis.matrix(1:2)[1] FALSEクラスも定義されています。
class(m)[1] "matrix" "array" 2.2 行名と列名の変更・削除
行名と列名は変更できます。
(rownames(m) <- c("r1", "r2"))[1] "r1" "r2"(colnames(m) <- c("c1", "c2"))[1] "c1" "c2"NULLを代入したら削除できます。
3 行列の要素の参照と代入
行列の要素は行番号と列番号はもちろん、行名と列名でも参照可能です。
# 行列を作成
# lettersはアルファベットが並んだビルトイン配列、toupperは大文字にする関数
(m <- matrix(1:9, 3, dimnames = list(letters[1:3], toupper(letters[1:3])))) A B C
a 1 4 7
b 2 5 8
c 3 6 9m[2, 3] # 2行3列を参照[1] 8m["b", "C"] # 行名と列名でも参照可[1] 8m[1, ] # 1行目すべてをベクトルとして参照A B C
1 4 7 m["a", ] # 行名でも参照可能A B C
1 4 7 m[, 3] # 3列名すべてをベクトルとして参照 a b c
7 8 9 m[, "C"] # 列名でも参照可能a b c
7 8 9 m[5] # ベクトルとして5番目の要素を参照[1] 5m[c(1,3), c(1,3)] # 1行目と3行目、1列目と3列名に含まれる小行列 A C
a 1 7
c 3 9c(m) # 行列の成分をベクトルにする[1] 1 2 3 4 5 6 7 8 9参照と同様に代入することができます。
m[1, 1] <- 8
m[1, c(2, 3)] <- c(3, 4)
m[, 2] <- c(3, 9, 6)
m[3, ] <- c(5, 6, 1)
m["b", "C"] <- 7
m A B C
a 8 3 4
b 2 9 7
c 5 6 1m[c(1,3), c(1,3)] <- 0
m A B C
a 0 3 0
b 2 9 7
c 0 6 0# ベクトルと見なして9番目の要素にNAを代入
m[9] <- NA
m A B C
a 0 3 0
b 2 9 7
c 0 6 NANAを入れましたが、Infなども入ります。
diagを使うと対角成分に代入もできます。
diag(m) <- c(1.1, 2.2, 3.3)
m A B C
a 1.1 3.0 0.0
b 2.0 2.2 7.0
c 0.0 6.0 3.33.1 行と列の削除
ベクトルと同様に行番号や列番号をマイナスにすると、その行や列を削除した行列になります。
m[-2, ] # 2行目を除いた小行列 A B C
a 1.1 3 0.0
c 0.0 6 3.3m[, -2] # 2列目を除いた小行列 A C
a 1.1 0.0
b 2.0 7.0
c 0.0 3.3m[-2, -2] # 2行目と2列目を除いた小行列 A C
a 1.1 0.0
c 0.0 3.3残す行と列のみを選択するのも方法です。
m[c(1, 3), ] # 1行目と3行目を残した小行列 A B C
a 1.1 3 0.0
c 0.0 6 3.3m[, c(1, 3)] # 1列目と3列目を残した小行列 A C
a 1.1 0.0
b 2.0 7.0
c 0.0 3.3m[c(1, 3), c(1, 3)] # 1, 3行目と1, 3列目を残した小行列 A C
a 1.1 0.0
c 0.0 3.34 行列の演算
統計学の教科書で使うような演算ができます。
# 同じ大きさの行列の四則演算
matrix(1:9, 3, 3) + diag(3) [,1] [,2] [,3]
[1,] 2 4 7
[2,] 2 6 8
[3,] 3 6 10# 転置
(m <- matrix(1:4, 2)) [,1] [,2]
[1,] 1 3
[2,] 2 4t(m) [,1] [,2]
[1,] 1 2
[2,] 3 4# 行列積
matrix(1:6, 2) %*% matrix(-1:-6, ncol = 2) [,1] [,2]
[1,] -22 -49
[2,] -28 -64# 逆行列
solve(m) [,1] [,2]
[1,] -2 1.5
[2,] 1 -0.5# クロネッカー積
diag(2) %x% matrix(1:4, 2, 2) [,1] [,2] [,3] [,4]
[1,] 1 3 0 0
[2,] 2 4 0 0
[3,] 0 0 1 3
[4,] 0 0 2 4# 外積
matrix(1:2, 2, 1) %o% matrix(-1:-2, 2, 1), , 1, 1
[,1]
[1,] -1
[2,] -2
, , 2, 1
[,1]
[1,] -2
[2,] -4なお、solveはX %*% beta = yのとき、solve(X, y)としたらbetaを出す連立方程式を計算する関数で、yのデフォルトが単位行列になっています。
上三角を操作したい場合は、upper.triを使います。
4.1 行列とベクトルの演算
行列をベクトルとして扱った後、行列に戻されます。
# 行列の各要素への演算
matrix(1:9, 3, 3) %/% 5 # 5で割って端数切捨て
# ベクトルを足す
# matrix(1:9 + c(1, 3, 2), 3)と同じ結果
matrix(1:9, 3) + c(1, 3, 2)4.2 下三角と上三角
まとめて操作するための関数があります。真偽値を戻し、diagで対角成分を含めるか否かを指定できます。
X <- matrix(1:25, 5, 5)
X[!lower.tri(X, diag = TRUE)] <- 0
X [,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 2 7 0 0 0
[3,] 3 8 13 0 0
[4,] 4 9 14 19 0
[5,] 5 10 15 20 254.3 行列式
行列式を求める関数が用意されています。
(X <- matrix(round(runif(25, min=0.5, max=9.5)), 5)) [,1] [,2] [,3] [,4] [,5]
[1,] 8 8 5 4 2
[2,] 4 3 3 1 6
[3,] 9 5 2 2 5
[4,] 6 4 1 7 8
[5,] 6 7 9 2 1det(X)[1] -2787もちろん正方行列にしか使えません。
4.4 固有値と固有ベクトル
固有値と固有ベクトルの計算をする関数も用意されています。
(e <- eigen(X))eigen() decomposition
$values
[1] 23.4252076+0.00000i -6.8886359+0.00000i 3.6843904+0.00000i
[4] 0.3895189+2.12977i 0.3895189-2.12977i
$vectors
[,1] [,2] [,3] [,4]
[1,] 0.4968143+0i 0.3570585+0i 0.17618808+0i -0.44103229-0.2316812i
[2,] 0.3191065+0i -0.3621329+0i -0.40922203+0i 0.65900103+0.0000000i
[3,] 0.4360642+0i -0.4297990+0i -0.06907061+0i -0.28667312+0.2230519i
[4,] 0.5062277+0i -0.3865755+0i 0.85035627+0i -0.07909817-0.2304676i
[5,] 0.4526899+0i 0.6381232+0i -0.27132740+0i 0.16382283+0.3152595i
[,5]
[1,] -0.44103229+0.2316812i
[2,] 0.65900103+0.0000000i
[3,] -0.28667312-0.2230519i
[4,] -0.07909817+0.2304676i
[5,] 0.16382283-0.3152595iprod(e$values) # 固有値をすべてかければ行列式になる[1] -2787+0i5 行列の中から指定した要素の位置を特定する
行列の場合はarr.ind=TRUEオプションをつけておくと、whichで行番号と列番号をそれぞれ得られます。
(w <- matrix((1:12)[order(runif(12))], 4, 3)) [,1] [,2] [,3]
[1,] 5 1 7
[2,] 3 12 8
[3,] 9 2 10
[4,] 6 11 4(i <- which(w == 6, arr.ind = TRUE)) row col
[1,] 4 1w[i][1] 66
配列（array）
行列はベクトルを2次元配列と見なすものでしたが、3次元以上の配列と見なす事もできます。
(a <- array(1:24, c(2, 3, 4))), , 1
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
, , 2
[,1] [,2] [,3]
[1,] 7 9 11
[2,] 8 10 12
, , 3
[,1] [,2] [,3]
[1,] 13 15 17
[2,] 14 16 18
, , 4
[,1] [,2] [,3]
[1,] 19 21 23
[2,] 20 22 24array(1:24, c(2, 3, 4), dimnames = list(c("a", "b"), c("A", "B", "C"), c("alpha", "beta", "gamma", "delta"))), , alpha
A B C
a 1 3 5
b 2 4 6
, , beta
A B C
a 7 9 11
b 8 10 12
, , gamma
A B C
a 13 15 17
b 14 16 18
, , delta
A B C
a 19 21 23
b 20 22 24行列と同様に参照と代入と演算ができます。%*%や%x%や%o%も定義されています。また、多次元配列のうち2次元を取り出すと、行列になります。