R における「配列 (array)」は、ベクトルを多次元に拡張したものである。2次元のものは特に「行列 (matrix)」と呼ぶ。逆に言うと、R では3次元以上のものは行列とは呼ばない。
配列
配列の実体は、次元ベクトル dim を属性に持つベクトルである。次元ベクトルとは、正整数の数値ベクトルでできている。次元ベクトルの要素数が配列の次元数であり、次元ベクトルの各要素の値は、各次元の要素数 (添字の上限値) を示す。
生成
例えば 2×3 の2次元配列を作るには、以下のような方法がある。
配列を生成する関数 array() を用いた方法
> a <- array(0, dim=c(2, 3))
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
array() 関数を用いると、データベクトルと次元ベクトルから配列が生成できる。データベクトルの要素数が次元ベクトルから計算される総要素数に足りないときは、例によってデータベクトルの要素が先頭から再利用される。データベクトルの方が大きいときは必要なところまでが使われる。分かりやすい例としては、array(0, c(3, 4, 2)) とすれば、全要素が 0 の配列になる (要素数1のベクトルが必要なだけ再利用されているということ)。
2×3 = 6 要素のベクトルに次元ベクトルを付加する方法
> a <- rep(0, 6)
> dim(a) <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
dim 属性を dim() ではなく attr() で直接操作する方法
出来るというだけで、当然だがおすすめしない。
> a <- rep(0, 6)
> attr(a, "dim") <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
添字の順番
ベクトルに次元ベクトルを付加して配列とした場合に、添字が割り振られる順番を確認してみる。
> a <- 1:6
> dim(a) <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
a[1,1] a[2,1] a[1,2] a[2,2] a[1,3] a[2,3] の順に割り当てられていることが分かる。大雑把に言えば、最初の添字が先に変わる順番。これは FORTRAN と同じで、C, Java とは異なる。
添字とサブセクション
配列の中身を参照するにあたっては、単純な正整数値の指定だけではなく、添字ベクトルが有効である。次元数分添字ベクトルを「,」区切りで並べる。空の添字ベクトルが与えられた次元は、無条件で全てを対象とする意味になる (TRUE を指定したのと同じ)。
> a[1,]
[1] 1 3 5
> a[2,]
[1] 2 4 6
> a[,1]
[1] 1 2
> a[,2]
[1] 3 4
> a[,3]
[1] 5 6
> a[,]
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
a[] や a[,] は配列全体の意味になり、つまり a と同じである。
添字に単一かつ配列ではないベクトルが指定された場合は、次元ベクトルが無視されデータベクトルに対する参照となる。配列を指定した場合は、以下に述べる。
添字配列
添字として、添字ベクトルの代わりに「添字配列」を用いることができる。不規則に配列中の位置を指定して値を取得したり代入したりできる。
例えば、4×5 の2次元数値配列 X の中の X[1,3]、X[2,2]、X[3,1] の値を 0 に置き換えたいとする。X は以下のようなもの。
> X <- array(1:20, dim=c(4, 5))
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 5 9 13 17
[2,] 2 6 10 14 18
[3,] 3 7 11 15 19
[4,] 4 8 12 16 20
操作するにあたって、[1, 3], [2, 2], [3, 1] を意味する添字配列 i を作る。
> i <- array(c(c(1, 2, 3), c(3, 2, 1)), dim=c(3, 2))
> i
[,1] [,2]
[1,] 1 3
[2,] 2 2
[3,] 3 1
X[i] を見てみる。
> X[i] [1] 9 6 3
X[i] を 0 で置換する。
> X[i] <- 0
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 5 0 13 17
[2,] 2 0 10 14 18
[3,] 0 7 11 15 19
[4,] 4 8 12 16 20
より複雑な例。blocks, varieties の2つの因子を持つブロック実験計画用の計画行列を作成する。実験はn回の反復がある。blocks は3水準、varieties は2水準、実験は12回とする。
> blocks <- factor(1:3) > varieties <- factor(1:2) > n <- 12
blocks, varieties のそれぞれについて、実験回数×水準数の行列を作り、Xb, Xv とする。
> Xb <- matrix(0, n, nlevels(blocks)) > Xv <- matrix(0, n, nlevels(varieties))
blocks, varieties のそれぞれについて、実験に水準を順番に割り当てた配列 ib, iv を作る。
> ib <- cbind(1:n, blocks)
> iv <- cbind(1:n, varieties)
> ib
blocks
[1,] 1 1
[2,] 2 2
[3,] 3 3
[4,] 4 1
[5,] 5 2
[6,] 6 3
[7,] 7 1
[8,] 8 2
[9,] 9 3
[10,] 10 1
[11,] 11 2
[12,] 12 3
> iv
varieties
[1,] 1 1
[2,] 2 2
[3,] 3 1
[4,] 4 2
[5,] 5 1
[6,] 6 2
[7,] 7 1
[8,] 8 2
[9,] 9 1
[10,] 10 2
[11,] 11 1
[12,] 12 2
ib, iv は、1列目に実験の回号、2列目に水準の入った2次元配列である。
ib, iv を添字配列として、Xb, Xv の該当箇所に 1 を入れる。
> Xb[ib] <- 1
> Xv[iv] <- 1
> Xb
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
[4,] 1 0 0
[5,] 0 1 0
[6,] 0 0 1
[7,] 1 0 0
[8,] 0 1 0
[9,] 0 0 1
[10,] 1 0 0
[11,] 0 1 0
[12,] 0 0 1
> Xv
[,1] [,2]
[1,] 1 0
[2,] 0 1
[3,] 1 0
[4,] 0 1
[5,] 1 0
[6,] 0 1
[7,] 1 0
[8,] 0 1
[9,] 1 0
[10,] 0 1
[11,] 1 0
[12,] 0 1
これで、Xb, Xv は各実験で用いる水準をフラグ 1 で示した配列となった。
> X <- cbind(Xb, Xv)
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 1 0
[2,] 0 1 0 0 1
[3,] 0 0 1 1 0
[4,] 1 0 0 0 1
[5,] 0 1 0 1 0
[6,] 0 0 1 0 1
[7,] 1 0 0 1 0
[8,] 0 1 0 0 1
[9,] 0 0 1 1 0
[10,] 1 0 0 0 1
[11,] 0 1 0 1 0
[12,] 0 0 1 0 1
Xb, Xv を併せれば、求める計画行列になる。因子に名前を付けていないため分かりづらいが、[,1], [,2], [,3] が blocks、[,4], [,5] が varieties に対応している。
生起行列は Xb, Xv を用いれば、以下のように求まる。
> N <- crossprod(Xb, Xv)
> N
[,1] [,2]
[1,] 2 2
[2,] 2 2
[3,] 2 2
table() を用いれば、ib, iv から直接求めることもできる。
> N <- table(ib[,2], iv[,2])
> N
1 2
1 2 2
2 2 2
3 2 2
配列の演算
配列は数値演算式で用いることができ、データベクトルの要素毎に演算を行った結果の配列となる。基本的には被演算項も同じ次元ベクトルを持つ必要がある。
2つの配列の外積
配列に対する重要な演算として「外積 (outer product)」がある。a と b の外積は、次元ベクトルとして a と b の次元ベクトルを連結したものを持ち、データベクトルには a と b の各要素全ての組み合わせの積が格納される。外積は演算子「%o%」または関数 outer() で演算子「*」を用いることで表現する。
> a <- array(1:2, c(1, 2)) > b <- array(1:12, c(3, 4)) > ab <- a %o% b > ab <- outer(a, b, "*") > dim(a) [1] 1 2 > dim(b) [1] 3 4 > dim(ab) [1] 1 2 3 4
例として、各要素が 0～9 の整数値を取る2×2行列 $\left(
\begin{array}{cc}
a & b \\
c & d
\end{array}
\right)$ の行列式 $ad - bc$ の値について、各要素の値が独立に一様に選ばれた場合の行列式の値の分布を求める。
> d <- outer(0:9, 0:9)
> d
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 2 3 4 5 6 7 8 9
[3,] 0 2 4 6 8 10 12 14 16 18
[4,] 0 3 6 9 12 15 18 21 24 27
[5,] 0 4 8 12 16 20 24 28 32 36
[6,] 0 5 10 15 20 25 30 35 40 45
[7,] 0 6 12 18 24 30 36 42 48 54
[8,] 0 7 14 21 28 35 42 49 56 63
[9,] 0 8 16 24 32 40 48 56 64 72
[10,] 0 9 18 27 36 45 54 63 72 81
ベクトル 0:9 と 0:9 の outer() をとることで、0*0～9*9 の全ての値が含まれる配列ができる。これが $ad$ および $bc$ の値の候補になる。
> fr <- table(outer(d, d, "-")) > fr -81 -80 -79 -78 -77 -76 -75 -74 -73 -72 -71 -70 -69 -68 -67 -66 -65 -64 -63 -62 -61 -60 -59 -58 -57 -56 -55 -54 -53 -52 -51 -50 19 1 2 2 3 2 4 2 4 41 4 4 8 6 6 10 7 27 49 8 8 17 8 12 18 53 13 60 12 18 22 16 -49 -48 -47 -46 -45 -44 -43 -42 -41 -40 -39 -38 -37 -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 35 70 22 24 66 28 18 72 22 75 37 34 26 111 63 36 45 84 34 94 36 93 97 50 53 156 42 60 103 107 50 168 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 51 140 112 116 59 191 65 126 156 185 115 206 117 179 153 156 111 570 111 156 153 179 117 206 115 185 156 126 65 191 59 116 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 112 140 51 168 50 107 103 60 42 156 53 50 97 93 36 94 34 84 45 36 63 111 26 34 37 75 22 72 18 28 66 24 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 22 70 35 16 22 18 12 60 13 53 18 12 8 17 8 8 49 27 7 10 6 6 8 4 4 41 4 2 4 2 3 2 79 80 81 2 1 19
d と d の outer() を演算子「-」でとることで、$ad - bc$ としてありうる全ての値がとれる。table() で度数を求めている。上記の例では分かりにくいが、奇数行が出現した値で、偶数行はその上の行の値が出現した回数が表示されている。これを plot() でヒストグラムにするには、以下のようにする。
> plot(as.numeric(names(fr)), fr, type="h", xlab="Determinant", ylab="Frequency")
X 軸の値を as.numeric() で数値に強制変換しているのは、出現した数値が fr 内では names 属性の値として文字列になってしまっているため。仮にループで解こうとすると、以下のようになる。
> fr <- integer(); for (a in 0:9) for (b in 0:9) for (c in 0:9) for (d in 0:9) { det <- as.character(a*d - b*c); fr[det] <- ifelse(is.na(fr[det]), 0, fr[det]) + 1 }
配列の一般化転置
aperm(a, perm) は置換ベクトル perm を用いて a の次元を入れ替える。perm は 1～dim(a) の置換になっている必要がある。新しい配列の j 番目の次元に perm[j] を持ってくる。行列の転置の一般化である。つまり、行列 A の転置は aperm(A, c(2, 1)) で得られる。
> A <- matrix(1:6, c(2, 3))
> A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> aperm(A, c(2, 1))
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
行列
「行列 (matrix)」とは R においては2次元の配列のことを指す。R ではその他の次元のものは行列とは呼ばない。行列には行列専用の便利な演算、関数が定義されている。
転置、行数、列数
t() は行列の転置を行う。任意の次元の配列における一般化転置は aperm() を用いるが、行列においては t() で簡便に行うことができる。nrow(), ncol() はそれぞれ行数、列数を取得する関数である。これらは配列に対しても1次元目、2次元目の要素数を返す関数として機能する。
作成
行列は matrix() で作る。
> A <- matrix(c(c(1, 2), c(3, 4)), 2)
> B <- matrix(c(c(0, -1), c(-1, 0)), 2)
> A
[,1] [,2]
[1,] 1 3
[2,] 2 4
> B
[,1] [,2]
[1,] 0 -1
[2,] -1 0
積
演算子「%*%」は行列の積を求める演算子である。ベクトルも与えることができて適宜行ベクトルまたは列ベクトルとして解釈される。n×1行列、1×n行列も行ベクトルまたは列ベクトルとして扱われる。ベクトルが行ベクトルか列ベクトルかの解釈は、例えば x %*% x といった式では曖昧になる。t(x) %*% x と x %*% t(x) のどちらに変換すべきかは自明ではない。R では、より小さい行列、すなわち t(x) %*% x のスカラー値が結果となる。他方の値を求めるには x %*% t(x) の他に cbind(a) %*% a や a %*% rbind(a) のようにいずれかのオペランドを行列として明示しても良い。もちろん cbind(a) %*% rbind(a) でも良い。
> A * B
[,1] [,2]
[1,] 0 -3
[2,] -2 0
A * B は単に成分毎の積からなる行列であり、行列の積ではない。
> A %*% B
[,1] [,2]
[1,] -3 -1
[2,] -4 -2
%*% で行列の積が求まる。
> A %*% B[,1]
[,1]
[1,] -3
[2,] -4
> A %*% B[,2]
[,1]
[1,] -1
[2,] -2
> A[1,] %*% B
[,1] [,2]
[1,] -3 -1
> A[2,] %*% B
[,1] [,2]
[1,] -4 -2
%*% で行列とベクトルの積も計算できる。
クロス積
crossprod() 関数は行列のクロス積を計算する。crossprod(A, B) は t(A) %*% B と等しいが、より効率的に動作する。
> crossprod(A, B)
[,1] [,2]
[1,] -2 -1
[2,] -4 -3
> t(A) %*% B
[,1] [,2]
[1,] -2 -1
[2,] -4 -3
その他行列演算
diag() 関数は、ベクトルを渡すとそのベクトルの値を対角成分に持つ行列を返す。行列を渡すと対角成分を抜き出したベクトルを返す。Matlab と同様の挙動。正整数 k を引数に渡すと、k×kの単位行列を返す。ややこしい。
solve(A, b) 関数は線形方程式 $Ax = b$ を解いて x を求める。
> A <- matrix(c(3, 2, 4, 5), 2) > b <- c(11, 12) > solve(A, b) [1] 1 2
solve(A) は A の逆行列を求める。しかし、solve(A, b) は solve(A) %*% b よりも効率が良いので、二次形式を求めるときも x %*% solve(A) %*% x とするより x %*% solve(A, x) とした方が効率が良い。solve(A) を用いることは稀だ。
eigen(Sm) は対称行列 Sm の固有値と固有ベクトルを求める。固有値 values と固有ベクトル vectors のリストとして返ってくる。eigen(Sm)$values, eigen(Sm)$vectors などとすればそれぞれ固有値、固有ベクトルのみが取り出せる。$ は再帰的構造のオブジェクトから要素を取り出す演算子。
svd(M) は行列 M の特異値分解を行う。特異値分解とは、M を U %*% D %*% t(V) の形に分解することである。M が m 行 n 列の行列とすると、U は m 行 m 列のユニタリ行列、V は n 行 n 列のユニタリ行列、D は m 行 n 列で対角成分は非負、それ以外の成分は 0 の行列となる。svd(M) の結果は d, u, v の3つの要素を持つリストとして返される。ただし d は行列ではなく対角成分を並べたベクトルである (行列にしたければ diag(d))。M の行列式の絶対値は prod(svd(M)$d) で求まる。
正方行列のトレース (対角和) tr() という関数を定義してみる。
> tr <- function(X) ifelse(is.matrix(X), sum(diag(X)), X[1]) > tr(A) [1] 8 > tr(9.9) [1] 9.9 > tr(-9) [1] -9 > tr(c(1,2,3,4,4) + ) [1] 1 > tr(integer()) [1] NA
行列が指定された場合は、diag() で対角成分を取り出し、sum() で総計する。それ以外が指定されたときはベクトルとみなして先頭の要素を返す。結果としてスカラーなら指定された値がそのまま返る。3次元以上の配列が指定されたときも先頭の要素が返るが、それが合理的かどうかはわからない。
関数 lsfit(X, y) は計画行列 X と観測値のベクトル y に対して最小二乗法によるフィッティングを行う (Least squares fit)。ヘルプに載っている例は以下のもの。
> ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
> trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
> weight <- c(ctl, trt)
> lsD9 <- lsfit(x = unclass(gl(2,10)), y = weight)
> ls.print(lsD9)
Residual Standard Error=0.6964
R-Square=0.0731
F-statistic (df=1, 18)=1.4191
p-value=0.249
Estimate Std.Err t-value Pr(>|t|)
Intercept 5.403 0.4924 10.9723 0.000
X -0.371 0.3114 -1.1913 0.249
lsfit(X, y) は内部で QR 分解を行っている。
> X <- cbind(1, unclass(gl(2, 10))) > Xplus <- qr(X)
この例は、計画行列の QR 分解を行い、Xplus に付値している。
qr(X) で X の QR 分解を得る。Xplus はクラス qr で list モードのオブジェクト。qr, rank, qraux, pivot という要素を持つ。qr(X)$qr は X と同じ次元の行列で、上三角部分が R になっており、残りの要素と qr(X)$qraux が Q の情報をコンパクトな形式で持っている。qr.X(X) で元の行列が取り出せ、qr.Q(X), qr.R(X) でそれぞれ分解された直交行列と上三角行列が取り出せる。当然、qr.Q(X) %*% qr.R(X) == qr.X(X) == X。
> b <- qr.coef(Xplus, y)
> fit <- qr.fitted(Xplus, y)
> res <- qr.resid(Xplus, y)
> b
X
5.403 -0.371
> fit
[1] 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661
> res
[1] -0.862 0.548 0.148 1.078 -0.532 -0.422 0.138 -0.502 0.298 0.108 0.149 -0.491 -0.251 -1.071 1.209 -0.831 1.369
[18] 0.229 -0.341 0.029
b, fit, res はそれぞれ、y の X の範囲への射影の係数ベクトル (coefficients)、直交射影、直交補空間への射影 (residuals)。X はフルランク (全ての行または列が線形独立) でなくても良く、冗長なものは取り除かれる。
cbind() と rbind()
ベクトルや行列を合併して新しい行列を作る関数として、cbind(), rbind() がある。与えられたベクトルや行列を順番に cbind() では列方向に (右に)、rbind() では行方向に (下に) 並べる。
cbind() の引数に行列を指定する場合は、指定した全ての行列の行数が等しい必要がある。ベクトルの要素数は任意でよく、足りなければ先頭から再利用される。行列が指定されていない場合は要素数が最大のベクトルの要素数を行数とした行列が作られる。短いベクトルの要素の再利用もその要素数まで行われる。ここでの再利用は、ターゲットの行数が行数の少ない行列の行数や要素数の少ないベクトルの要素数の正整数倍でなければならない。そうでない場合はエラーとなる。
rbind() も行と列の位置付けが逆なことを除いて cbind() と同様。
配列に c() を適用すると……
c() も配列に対して適用することができる。しかし、cbind(), rbind() が dim 属性を意識しているのに対し、c() は無視している。つまり、データベクトルそのままの素のベクトルとみなして処理する。配列 X に対して c(X) とするだけでもデータベクトルが取り出される。
因子から度数分布表を作る
table() 関数は因子を用いて度数分布表を作成する関数である。4章の例を用いると、以下のようにできる。
> table(statef) statef act nsw nt qld sa tas vic wa 2 6 2 5 4 2 5 4
この例では、tapply(statef, statef, length) と同じことが行われている。
2つの因子を扱うこともできる。
> table(statef, sex)
sex
statef f m
act 2 0
nsw 2 4
nt 2 0
qld 1 4
sa 1 3
tas 1 1
vic 4 1
wa 2 2
新たに収入の値で区分する因子 incomef を定義する。
> incomef <- factor(cut(incomes, breaks = 35 + 10 * (0:7)))
cut(x, breaks) は x に与えられたベクトルを breaks に与えられたベクトルの値を境界にして分解する関数。これを使って3つの因子を指定した場合は以下のような表示になる。
> table(statef, sex, incomef)
, , incomef = (35,45]
sex
statef f m
act 1 0
nsw 0 1
nt 0 0
qld 0 1
sa 0 0
tas 0 0
vic 1 0
wa 0 0
, , incomef = (45,55]
sex
statef f m
act 1 0
nsw 1 0
nt 1 0
qld 0 1
sa 1 1
tas 0 0
vic 1 0
wa 2 1
, , incomef = (55,65]
sex
statef f m
act 0 0
nsw 1 2
nt 1 0
qld 1 2
sa 0 2
tas 1 1
vic 1 1
wa 0 1
, , incomef = (65,75]
sex
statef f m
act 0 0
nsw 0 1
nt 0 0
qld 0 0
sa 0 0
tas 0 0
vic 1 0
wa 0 0