Rで鶴亀算を解いてみる

工学や経済学の教育を受けてない人(例：筆者)は線形代数を使って現実の問題を解く機会に恵まれない。若くない人(例：筆者)は教科書を開いてみる前にまず自分が何を学ぶのかを知りたがる。ともかく誰でも知っている2行2列の最単純のケースでやってみよう。

Rで鶴亀算を解いてみる。

鶴と亀を混合し、頓着せずに数えると頭が3個、足が10本。
関数cは要素からベクトルの生成。関数rbindはベクトルを行ベクトルとして行列を生成。

> b = c(3,10) ; b
[1]  3 10
> head = c(1,1)
> feet = c(2,4)
> m = rbind(head,feet); m
     [,1] [,2]
head    1    1
feet    2    4
> solve(m, b)
[1] 1 2

鶴が1羽、亀が2匹。内部的には消去法を使っているのだろうか？
(追記)R 2.14.0のsolve.Rを読むと通常はLAPACKのLA_GESVを呼び出す。LAPACKでの実装はLU分解によるGauss消去法。ただしRのsolveからはpivot選択の結果は取れない

わざわざ反復法を使って解いてみる。まずJacobi法。

係数行列mを対角要素dとそれ以外rに分け、dの逆行列を求めて式にぶちこみ、回す。2x2行列では全くの無駄。行列やベクトルの掛け算を %*% と書くのが戸惑う。

> b = c(3,10) ; b
[1]  3 10
> d_1 = c(1,0)
> d_2 = c(0,4)
> r_1 = c(0,1)
> r_2 = c(2,0)
> d = rbind(d_1,d_2); d
d
    [,1] [,2]
d_1    1    0
d_2    0    4
> r = rbind(r_1,r_2); r
    [,1] [,2]
r_1    0    1
r_2    2    0
> d_inv = solve(d); d_inv
     d_1  d_2
[1,]   1 0.00
[2,]   0 0.25
> i = c(0,0)
> x_ = i
> x = d_inv %*% (b - r %*% x_) ; x_ = x ; x
     [,1]
[1,]  3.0
[2,]  2.5
> x = d_inv %*% (b - r %*% x_) ; x_ = x ; x
     [,1]
[1,]  0.5
[2,]  1.0

[1,] 2.00
[2,] 2.25

[1,] 0.75
[2,] 1.50

[1,] 1.500
[2,] 2.125

[1,] 0.875
[2,] 1.750

[1,] 1.2500
[2,] 2.0625

[1,] 0.9375
[2,] 1.8750

[1,] 1.12500
[2,] 2.03125

[1,] 0.96875
[2,] 1.93750

[1,] 1.062500
[2,] 2.015625

[1,] 0.984375
[2,] 1.968750

[1,] 1.031250
[2,] 2.007812

[1,] 0.9921875
[2,] 1.9843750

[1,] 1.015625
[2,] 2.003906

[1,] 0.9960938
[2,] 1.9921875

[1,] 1.007812
[2,] 2.001953

[1,] 0.9980469
[2,] 1.9960938

[1,] 1.003906
[2,] 2.000977

[1,] 0.9990234
[2,] 1.9980469

初期値(0,0)からの反復計算は、最初の結果(1,2)に収束していっている。

今度はGauss-Seidel法。

計算を要素毎にやって、計算の済んだ要素を別の要素の計算ですぐ使う点がJacobi法と違うところで、だから式としては行列の掛け算ではなくて、それぞれの要素の式になる。

Jacobiの時と明らかに収束の仕方が違うのだが、自分にはまだこの原因がわかってない。

(追記)逐次代入で一個前の値は本当は要らないし = は <- と書くのがRとしては正統派らしい

> b = c(3,10) ; b
[1]  3 10
> x_ = 0 ; y_ = 0
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 3 1
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 2.0 1.5
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.50 1.75
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.250 1.875
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.1250 1.9375
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.06250 1.96875
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.031250 1.984375
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.015625 1.992188
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.007812 1.996094
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.003906 1.998047
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.001953 1.999023
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000977 1.999512
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000488 1.999756
> x = 3 - y_ ; x_ = x ; y = .25 * (10 - 2 * x_) ; y_ = y ; c(x,y)
[1] 1.000244 1.999878

行列式を使う解法(Cramer's rule)

係数行列mと頭足ベクトルbは前と同じ。bでmの列を置き換え行列式を求め、係数行列の行列式det(m)で割る。

b は c(3,10) としても rbind(3,10) としても結果は変わらない。Rのベクトル(前者)は縦横の区別がなく、文脈によって適切に転置(キャスト)されるが、後者は2行1列の行列ということになり、これは自動で転置できない。なお当然ながら t(cbind(3,10)) としても同じ結果になる。

> b = c(3,10) ; b
[1]  3 10
> m = rbind(c(1,1),c(2,4)) ; m
     [,1] [,2]
[1,]    1    1
[2,]    2    4
> det(m)
[1] 2
> m_1 = cbind(b,c(1,4)) ; m_1
      b  
[1,]  3 1
[2,] 10 4
> det(m_1)
[1] 2
> m_2 = cbind(c(1,2), b) ; m_2
        b
[1,] 1  3
[2,] 2 10
> det(m_2)
[1] 4
> c(det(m_1) / det(m), det(m_2) / det(m))
[1] 1 2

自分のためのおまけ

• ls(); で定義済みオブジェクトの一覧(実体はobjects関数と同じっぽい)
• rm(object); で特定のオブジェクト削除
• 某設定ファイルのパーザ書こう。。
• 固有値とか

> m
     [,1] [,2]
[1,]    1    1
[2,]    2    4
> eigen(m)
$values
[1] 4.5615528 0.4384472

$vectors
           [,1]       [,2]
[1,] -0.2703230 -0.8719282
[2,] -0.9627697  0.4896337

> a = eigen(m)$vectors[,1]
> m %*% a
          [,1]
[1,] -1.233093
[2,] -4.391725
> eigen(m)$values[1] * eigen(m)$vectors[,1]
[1] -1.233093 -4.391725
> dim(m)
[1] 2 2
> rank(m)
[1] 1.5 3.0 1.5 4.0
> rownames(m) = c("うえ", "した")
> m
     [,1] [,2]
うえ    1    1
した    2    4
> colnames(m) = c("みぎ", "ひだり")
> m
     みぎ ひだり
うえ    1      1
した    2      4
> "これ"
[1] "これ"
> c("これ")
[1] "これ"
> m
     みぎ ひだり
うえ    1      1
した    2      4
> m[1,1] = 5
> m
     みぎ ひだり
うえ    5      1
した    2      4
> m[1,1]
[1] 5
> m[2,2]
[1] 4
> m[1]
[1] 5
> m[4]
[1] 4
> length(m)
[1] 4
> t <- matrix(c(1,5,3,7),ncol=2,byrow=TRUE)
> t
     [,1] [,2]
[1,]    1    5
[2,]    3    7
> class(t)
[1] "matrix"
> t <- as.table(t)
> class(t)
[1] "table"
> t
    A   B
A   1   5
B   3   7
> attributes(t)
$dim
[1] 2 2

$dimnames
$dimnames[[1]]
[1] "A" "B"

$dimnames[[2]]
[1] "A" "B"


$class
[1] "table"