2008年8月11日月曜日

An Introduction to R を読む (3) - 配列と行列

An Introduction to R を読む (2) – オブジェクトとカテゴリー (factor) の続き。

 

配列

ベクトルに dim 属性を与えると、配列 (array) になる。 dim 属性はベクトルで与えられ、次元ベクトル (dimention vector) と呼ばれる。この次元ベクトルの大きさ k に応じて、配列は k 次元であると言う。

ベクトル → dim を設定することにより → 配列

試してみる。

> a <- 1:24
> a
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> dim(a) <- c (4,6)
> a
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    5    9   13   17   21
[2,]    2    6   10   14   18   22
[3,]    3    7   11   15   19   23
[4,]    4    8   12   16   20   24

> dim(a) <- c(3,4,2)
> a
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

, , 2

     [,1] [,2] [,3] [,4]
[1,]   13   16   19   22
[2,]   14   17   20   23
[3,]   15   18   21   24

> dim(a) <- c(24)
> a
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

次元ベクトルが 2つの場合、行列 (matrix) 。 1 つ場合は、ベクトルと同じように扱える。 (多少違いがあるみたい。) 上記の例を見るとわかるように、ベクトルの先頭から、列方向に配置されていく。

 

要素の取得

配列の中の要素を取得する。先ほどの例の中で、次元ベクトルが c(3,4,2) の配列を対象とする。

> a[,,1]
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> a[,2,1]
[1] 4 5 6
> a[3,2,1]
[1] 6

> # [1,]の行を取得する
> a[1,,]
     [,1] [,2]
[1,]    1   13
[2,]    4   16
[3,]    7   19
[4,]   10   22
> a[1,2,]
[1]  4 16
> a[1,2,3]

ちなみに a[,,] は全ての要素を取得することになる。インデックスが空になっているところは、その要素全てを対象にしていることを表わす。

 

要素の設定

配列のインデックスに配列を指定して、要素を取得、値の設定をすることができる。

> a
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    4    7   10   13   16   19   22
[2,]    2    5    8   11   14   17   20   23
[3,]    3    6    9   12   15   18   21   24

> a[array(c(c(1,2,3,2,1,2,3,2),1:8), dim=c(8,2))] <- 0
> a
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    0    4    7   10    0   16   19   22
[2,]    2    0    8    0   14    0   20    0
[3,]    3    6    0   12   15   18    0   24

ここで指定したいインデックスについて考えるとき、上記の場合であれば、まず行方向についてインデックスを見る。そして、次にそれに対応させるように列方向のインデックスについて考える。最後に、それらを組み合わせてインデックスを指定できる行列へとなるように次元の値を設定する。

 

array(ベクトル, ディメンション) による配列の作成

array() を利用して配列を生成することができる。ただし、 対象の要素が dim で指定した要素数に満たないときは、要素を満たすように要素の先頭から繰り返される。

> array(1:20, dim=c(6,4))
     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    2    8   14   20
[3,]    3    9   15    1
[4,]    4   10   16    2
[5,]    5   11   17    3
[6,]    6   12   18    4

 

outer product

outer product 。 全てのベクトルの要素を掛け合わせる。(R: Outer Product of Arrays)

> c(1,2,3) %o% c(4,5,6)
     [,1] [,2] [,3]
[1,]    4    5    6
[2,]    8   10   12
[3,]   12   15   18

> outer(c(1,2,3),c(4,5,6), "*")
     [,1] [,2] [,3]
[1,]    4    5    6
[2,]    8   10   12
[3,]   12   15   18

outer() には自分で定義した関数を渡すことができる。ただし、その関数は引数を 2 つ取ること。例えば、適当に関数を作って、それを outer() に渡してみる。 (rel. ループと再帰)

> f <- function(x,y) (x + y) / y

> a <- c(1,2,3)
> b <- c(4,5,6)

> outer(a,b,f)
     [,1] [,2]     [,3]
[1,] 1.25  1.2 1.166667
[2,] 1.50  1.4 1.333333
[3,] 1.75  1.6 1.500000

 

行列

行列の作成

matrix() によって行列を作成できる。

> matrix(1:9)
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7
 [8,]    8
 [9,]    9

># 行数と列数を指定
> matrix(1:9, nrow=3, ncol=3)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

># 行方向に行列を作成する
> matrix(1:9, nrow=3, ncol=3, byrow=TRUE)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9

 

A を上記の matrix() によって生成した行列とする。

> A * A
     [,1] [,2] [,3]
[1,]    1   16   49
[2,]    4   25   64
[3,]    9   36   81

># 行列の積
> A %*% A
     [,1] [,2] [,3]
[1,]   30   66  102
[2,]   36   81  126
[3,]   42   96  150


># quadratic form ってなんだろ? 二次形式???
> a <- c(1,2,3)
> a %*% A %*% a
     [,1]
[1,]  228

># クロス積
> crossprod(A,A)
     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

># 上記は以下と同じらしい。
> t(A) %*% A
     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

crossprod() とはクロス積 - Wikipedia のことか。

 

転置行列

転置行列 – Wikipedia によると、

mn 列の行列 A に対して A の (i, j) 要素と (j, i) 要素を入れ替えた nm 列の行列、つまり対角線で成分を折り返した行列のことである。

 

> A <- array(1:9,dim=c(3,3))
> A
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> aperm(A, c(2,1))
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9

aperm() の特殊なケースとして t() を用いて同様に、

> t(A)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9

 

対角要素、単位行列

対角行列 - Wikipedia によると

正方行列であって、その対角成分((i,i)-要素)以外が零であるような行列のことである。 (…)

対角行列の転置行列は同一である。そのため対角行列は対称行列でもある。

diag() 関数にベクトルを引数として与える。

> diag(c(1,2,3,4,5))
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    2    0    0    0
[3,]    0    0    3    0    0
[4,]    0    0    0    4    0
[5,]    0    0    0    0    5

># 行列を引数に与えると対角要素が得られる
> diag(matrix(1:9, nrow=3,ncol=3))
[1] 1 5 9

># 数値 k を引数に与えると k × k の単位行列が得られる
> diag(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

単位行列 - Wikipedia によると、

単位的環上で定義される同じ型の正方行列同士の、積演算における単位元のことである。単位行列の対角成分には 1 が並び、他は全て 0 となる:

 

行列と連立方程式

考え方は、線型方程式系 - Wikipedia を参考に。例えば、次の連立方程式を行列で解くとする。

2x + 3y = 15
 x +  y =  1

R では solve() を用いて次のようにする。

> A <- matrix(c(2,3,1,1),nrow=2,ncol=2,byrow=TRUE)
> A
     [,1] [,2]
[1,]    2    3
[2,]    1    1

> b <- matrix(c(15,1))
> b
     [,1]
[1,]   15
[2,]    1

> solve(A,b)
     [,1]
[1,]  -12
[2,]   13

逆行列を求めるには、上記の solve() を用いて、

> A <- matrix(1:4, nrow=2, ncol=2)
> A
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> solve(A)
     [,1] [,2]
[1,]   -2  1.5
[2,]    1 -0.5

 

固有値と固有ベクトル

はぁ~未だよくわがんねぇ (@_@;)... パタッ(o_ _)o~†

とりあえず、イメージだけでもしておくかぁ。

結論を一言で言うと,”ほとんど” の線形写像はベクトルの”引き伸ばし(倍率が1以下ならば縮小)”と考えることができ,その引き伸ばしの方向を決めているのが固有ベクトルで倍率が固有値です。ただし,この様子は実数の世界で完全に捉えることは不可能で,複素数の世界において可能となります。

ときわ台学/固有値/固有ベクトルの幾何学的意味より

行列を 1個固定して考えてみます. この行列は何かよくわからないんですが線形変換を表します. たいていのベクトルはこの線形変換によって変な方向を向いてしまうんですが, まれに方向が変わらず長さだけが変わるベクトルがあります. このように「長さだけが変わるベクトル」がこの線形変換 (ひいては行列) の固有ベクトルとなります. で, 長さの変化率が固有値.

固有値、固有ベクトル、対角化...何のため? - 教えて!gooより

R では、

> Sm <- matrix(c(3,1,1,3),nrow=2,ncol=2)
> Sm
     [,1] [,2]
[1,]    3    1
[2,]    1    3
> eigen(Sm)
$values
[1] 4 2

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

># 個々の値を取得するには...
> ev%val
Error: unexpected input in "ev%val"
> ev$val
[1] 4 2
> ev$vec
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

 

その他

特異値分解 : svd(M)。lsfit() : 最小二乗法? lm(), qr()。 パタッ(o_ _)o~†

 

cbind(), rbind()

cbind(), rbind() は、それぞれ列方向、行方向にベクトルや行列をくっつけて行列を作る。これにより、ベクトルを行列の行または列のように扱うことができる。

> cbind(1,1:5)
     [,1] [,2]
[1,]    1    1
[2,]    1    2
[3,]    1    3
[4,]    1    4
[5,]    1    5

> cbind(1:3,1:5)
     [,1] [,2]
[1,]    1    1
[2,]    2    2
[3,]    3    3
[4,]    1    4
[5,]    2    5
Warning message:
In cbind(1:3, 1:5) :
  number of rows of result is not a multiple of vector length (arg 1)

> rbind(1,1:5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    2    3    4    5

> rbind(1:3,1:5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    1    2
[2,]    1    2    3    4    5
Warning message:
In rbind(1:3, 1:5) :
  number of columns of result is not a multiple of vector length (arg 1)


># 行列の場合
> A
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> cbind(A,1:5)
     [,1] [,2] [,3]
[1,]    1    3    1
[2,]    2    4    2
Warning message:
In cbind(A, 1:5) :
  number of rows of result is not a multiple of vector length (arg 2)
> rbind(A,1:5)
     [,1] [,2]
[1,]    1    3
[2,]    2    4
[3,]    1    2
Warning message:
In rbind(A, 1:5) :
  number of columns of result is not a multiple of vector length (arg 2)

 

次元のクリア

配列の次元をクリアする方法。 as.vector() と c()。

> a <- c(1,2,3,4,5,6); dim(a) <- c(2,3)
> a
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> b <- as.vector(a)
> b
[1] 1 2 3 4 5 6

> c <- c(a)
> c
[1] 1 2 3 4 5 6

 

クロス集計表

以下のように、クラス (class)  と性別 (sex)  を表わすベクトルがあるとする。これからクロス集計表を作りたい。

> class
[1] "A" "A" "B" "A" "C" "C" "B" "C" "C"
> sex
[1] "M" "M" "M" "M" "F" "F" "F" "F" "F"

> table(class,sex)
     sex
class F M
    A 0 3
    B 1 1
    C 4 0

cut() によって数値を特定の範囲に分類することができる。例えば、身長のベクトルがあり、これを 150cm から 5cm ごとの間隔で分ける場合、

> height <- c(170,168,175,160,163,153,155,170,165)

> cut(height, breaks = 150+5*(0:5)) -> heightf
> heightf
[1] (165,170] (165,170] (170,175] (155,160] (160,165] (150,155] (150,155]
[8] (165,170] (160,165]
Levels: (150,155] (155,160] (160,165] (165,170] (170,175]

上記の変数を使い、男女ごとに身長のクロス集計表を得る。

> table(heightf, sex)
           sex
heightf     F M
  (150,155] 2 0
  (155,160] 0 1
  (160,165] 2 0
  (165,170] 1 2
  (170,175] 0 1