文系のための「自由度調整済み決定係数」

ここまでの話では、重回帰分析の仕組みとモデルの読み方を整理してきた。
要するに、重回帰係数について、連立方程式を解き、
重回帰係数の値の正負と値の大小によって変数間の関係を観察する。

そうそう。特異値分解による擬逆行列から連立方程式を解く場合には、
以下のような式によって、重回帰係数を算出することができた。
少々、クドいかもしれないが、とりあえず再掲しておく。



一般的に、擬逆行列の解法は変数pが対象数nよりも小さい場合と（n > p）、
変数pが対象数nよりも大きい場合（n < p）で異なっているが、
特異値分解を用いた場合には近似的に計算することができる。 

とにかく、この式によって導出されたベクトルが重回帰係数であり、
この重回帰係数を用いて、次のような重回帰モデルを導出することができる。



さて、基本的な話は、これで収まっているのであるが、 問題となるのは、
このモデルが本当に正しいのか？という問題である。
まずは、このモデルの適合の度合いというか、説明力を知りたい。

最も単純な方法は、残差から重回帰モデルの適合性を評価する決定係数を用いる方法。
 予測値のバラツキ（回帰の平方和）「Sr」、残差のバラツキ（残差平方和）「Se」、
そして、目的変数の総平方和 「St」から以下の式によって表すことができる。



この式を変形し、以下のようにしたものが決定係数と呼ばれるものであった。



詳細に関しては、残差と決定係数の話ですでに整理済みの話である。
決定係数は、モデルで説明できない残差の大きさを反映した指標である。

ところで、この指標は直感的に理解しやすいのであるが、少々、問題がある。
実は、変数の数が増えれば増えるほど決定係数の値は高くなってしまうのである。
そこで登場するのが、「自由度修正済み決定係数」と呼ばれるものである。



ここでは、目的変数の総平方和「St」に対して対象数「n-1」で基準化し、
Seは、Stと同様に対象数で基準化すると同時に、説明変数の数「p」で基準化をしている。
つまり、全体では一つの対象における一つの変数当たりの残差でモデルを評価する。

このような調整を行うことで、対象数と変数の数の影響が除去される。
とにかく、前回のデータを用いて、自由度調整済み決定係数とやらを計算してみる。
まず、前回と同様に、インターネット経由で直接データを読み込む。
今回の演習で使用するデータも前回と同じ。

library(RCurl) 

ライブラリを読み込めたら、
以下のコマンドを実行する。

# Google Spreadsheet からデータをダウンロードする
data <- getURL("https://docs.google.com/spreadsheet/pub?key=0AtOtIs5BjRVhdDdyTmgtNm9DcEVJb3VLbnh3ZFlNb2c&single=true&gid=0&range=A1%3AH295&output=csv")

# ダウンロードしたデータを読み込む
X <- as.data.frame(read.table(textConnection(data), header=TRUE, row.names=1, sep=","))

次に、重回帰分析を実行する。前回と同様に、「墳長」を目的変数とし、
その他の変数から重回帰モデルを立てる。

# 念の為に最初のデータの状態を保存しておく。
O <- X

# 実際の目的変数も保存しておく。
y <- O[, 1]

 # 目的変数を元の行列から分離する
y <- X[, 1]
X <- as.matrix(X[, -1])

# 切片の計算のために先頭行を「1」の行列を付け加える。
Intercept <- matrix(rep(1, nrow(X)), ncol=1)
X <- cbind(Intercept, X)

# 特異値分解を行う
X.svd <- svd(X)
X.svd.v <- X.svd$v
X.svd.d <- diag(X.svd$d)
X.svd.u <- X.svd$u

 # 重回帰係数を求める
X.lm <- ((X.svd.v %*% solve(X.svd.d)) %*% t(X.svd.u)) %*% y 

# 解りやすいように名前を付ける。
rownames(X.lm) <- c("Intercept",colnames( X[,2:ncol(X)]))
colnames(X.lm) <- "coefficients" 

# 最後に重回帰モデルから予測する。
# なお、以下の計算はベクトルの成分の掛け算。
y_hat <- X.lm[1] + X.lm[2]*X[,2] + X.lm[3]*X[,3] + X.lm[4]*X[,4] + X.lm[5]*X[,5] + X.lm[6]*X[,6] + X.lm[7]*X[,7]

# まずは、StとSeを求めて、決定係数を算出する。
Se <- sum((y-y_hat)^2)
St <- sum((y-mean(y))^2)

# 以下が標準的な決定係数。
1-(Se/St) 

# 次に、自由度調整済み決定係数を求める
p <- ncol(X)-1
n <- length(y)
1-(Se/(n-p-1))/(St/(n-1)) 

実行結果は以下の通り。決定係数と自由度調整済み決定係数はかなり高い。 

> # 以下が標準的な決定係数。
> 1-(Se/St)  
[1] 0.9926216
> 
> # 次に、自由度調整済み決定係数を求める
> p <- ncol(X)-1
> n <- length(y)
> 1-(Se/(n-p-1))/(St/(n-1))  
[1] 0.9924674

決定係数は、モデルの当てはまりの良さを表す指標であるので、
今回の結果の場合は、モデル全体では99%以上の説明力を持っていることが解る。
したがって、今回のモデルは、非常に当てはまりの良いモデルだと言える。

さて、残差の検討もしておく。これも、必ず確認しておく必要がある。

たしか、残差プロットは、まんべんなく散らばっていた方が良かった。

残差に何らかの傾向が現れているということは、重回帰モデルとして問題がある。

# 残差を標準化してバラツキ具合を観察する。
plot(scale(y-y_hat), ylab="Residuals")
abline(h=0)
abline(h=c(-2, 2), col="gray", lty=2)

# 外れ値を検出する
(r.lower <- which(scale(y-y_hat) < -2)) # −2よりも小さいもの
(r.upper <- which(scale(y-y_hat) > 2))  # +2よりも大きいもの
1-((length(r.upper)+length(r.lower))/n)

そして、実行結果は以下の通り。

# 外れ値を検出する
> (r.lower <- which(scale(y-y_hat) < -2)) # −2よりも小さいもの
[1] 91 98 137 243 258
> (r.upper <- which(scale(y-y_hat) > 2)) # +2よりも大きいもの
[1] 53 94 121 127 131 162 164 263 264 267
> 1-((length(r.upper)+length(r.lower))/n)
[1] 0.9489796 

一般的に重回帰分析における残差は正規分布に従うとされている。
そして、標準正規分布に従うのであれば、±2σ に95.4%が入ることになるので、
簡易な方法として、標準化した値の±2の外側を「外れ値」とする。 

今回の場合、全体の約94.9%が、±2σの範囲内に収まっていて、
また、散布図を見る限りは、何の傾向も見えてこない。
したがって、やはり、このモデルは非常に良い状況を表していることになる。 

自由度調整済み決定係数を用いることで、モデルの適合の度合いを知ることができ、
さらに、F検定という方法を用いることでモデルの適合性を評価するができる。
これに関しては、次回に整理することにする。

文系のための「データの観察」（１）

ケトレーの話では、「平均人」というものが出てきた。
平均人は、対象の個性を測るための架空の擬似人格のようなものであった。
そして、個性というのは、この平均人からの距離（偏差）であった。

データの「真ん中」を見極めることと、
個性の強さとも言える「バラツキ」を評価すること、
この2つの基本的操作は、データ解析において、最も重要なことである。
特に、データの個性を表す「バラツキ」はデータ解析において基本中の基本。

近年、文系分野の研究者が「見様見真似」で難しい分析を行なっているの見かける。
確かに、私自身、そうした時期があったし、そういった悪しき癖は残っている。
しかしながら、基本的な事だけでも、重要なことは分かる。

最初に考えるべきことは、データの「真ん中」である。「バラツキ」はその後。
データの「真ん中」が分からないと、「バラツキ」を計るための基準も定まらない。
本当は「バラツキ」について考えたいのだが、まずは「真ん中」の話から。

そもそも、データの「真ん中」とは何か？

多くの人は、「平均」という言葉を思い浮かべるのではないか？
なるほど、確かに、平均というのは良い考えである。
実際、ケトレーは、平均という概念が、データ全体を説明するものと考えた。

では、平均とは何か？魚の値段を例に考えてみる。
まぁ、魚である必要は無いのだが、私は魚介類が好きなのだ。

仮に、魚屋さんで400円のお刺身、200円の目刺し、600円の鰈の一夜干し、
この3つの品目が並んでいたとする。
適当だが、妥当な値段である。まぁ、良い。単なる実験データ。

ふむ。では、この3つの品目の値段の平均は？電卓を使っても良い。
もちろん、Rで計算しても良い。計算式は、自分で立てられるだろう。

答えは、400円。さすがに、この計算は大丈夫であろう。
この程度の計算ができないと困る。「文系」とかそういう問題では無い。

では、どような計算を行ったのか？数式を見ながら考えてみる。
おそらく、以下のような計算を考えたはずである。



あるいは、以下のように考えた人もいるはず。



2つ目の式は、品目3つを一つの単位（つまり全体を「3」）としていて、
したがって、値段全体を単位一つ分として考えた場合の、
単位一つ分の値段が理屈上どうなるかを考えている。どこかで聞いたような？

確か、「文系のための「逆行列」（1）」の話。

なるほど、1/3 を掛けるというのは、
対象数に合致するような単位に換算するために、
対象数の逆数をかけているのだった。

要するに、3つの品目に個性が無く、全て同じ値段だったら？というわけである。
したがって、400円という単位が3つで全体が説明できる。

ところで、上記の場合は、たった三品であったが、
品目が多くなると、式が長すぎて書くのも、見るのも面倒。
そこで、上記の式を抽象化し、より一般的なものに書き直してみる。



まず、品目の変数をxと置く。右下の数字と記号は何であったか？
これは「添字」と呼んだ。何番目の対象かを示していて、
不特定番目の対象を「i」番目、
最後番目の対象を「n」番目としている。
したがって、全体の個数（総数）は当然「n」個。

なるほど、これで、対象がいくら増えたとしても、
この長さの式で平均を書き表すことができる。

だが、やはりまだ長い。しかも、「・・・」という記号は美しくない。
ということで、自称「文系」が苦手とする「Σ（シグマ）」を使って表現する。



イコールで分けられた3つの式は同じ。
左端は、nで割るというイメージ。自称文系には、解りやすいが見難い。
真ん中は、nの逆数を掛けるというイメージ。理屈っぽい式。教科書はこの形式。
右端は、慣例的な表現で、「エックスバー」と読む。玄人向け。計算式を省略。
上に「-」が付いているものは、何かの平均を表す。

さて、ここで、シグマの見方を復習しておく。
Σという記号は、「足し算の繰り返し」を表す。
Σの記号の「下」についてあるのが、添字の記号とその開始番号。
Σの記号の「上」についてあるのが、添字の最後（総数）を表している。

一回目の足し算が終わる度に、iの値が1〜nまで一つずつ繰り上がる。
そして、i=nとなったときに、足し算の計算が終了する。
その結果を総数nで割ると考えるか、総数nの逆数で掛けると考えるかで、
式の表現は異なる。ただし、結果は同じ。

「エックスバー」による表現は、非常によく出てくる。
平均というのは、誰もが理解している（という前提）ので、
慣例的にこの表現を使う。説明無しに登場することも多い。
この意味が解らないと、教科書の内容をサッパリ理解できない。

以上で基本的な「平均」の話は終わり。

さて、ここからは、「R」を使って平均の計算をやってみる。
今度は、「魚の卸売価格」を例に考えてみる。
私は、魚介類が好きでなのである。

まずは、コンソールを立ち上げて、データ読み込みのコマンドを入力。実行はまだ。

データの読み込みは、read.table() 関数を使うのだった。

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",", header=TRUE)

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",", header=TRUE)

次に、今回の例で用いるデータをコピーした後、上のコマンドを実行する。

Oroshi_Kg, Taka, Naka, Yasu
Tai, 1336, 3150, 1217, 53
Chinu, 226, 630, 513, 105
Sawara, 212, 1575, 1259, 840
Yazu, 4834, 840, 315, 179
Suzuki, 618, 1575, 1146, 525
Akou, 21, 4725, 3129, 1575
Kochi, 28, 2100, 874, 210
Okoze, 28, 5250, 2635, 210
Ainame, 29, 3150, 621, 315
Hage, 1205, 1890, 383, 210
Konoshiro, 21, 2100, 1100, 105
Sayori, 12, 5250, 3916, 1260
Anago, 1637, 2625, 1721, 630
Mebaru, 330, 4200, 1680, 315
Tachiuo, 1211, 2310, 930, 105
Hamo, 1622, 3938, 592, 53
Karei, 41, 6300, 3178, 525
Hirame, 540, 2100, 1496, 210
Aji, 5149, 3150, 789, 210
Saba, 5413, 3675, 625, 53
Ika, 2414, 2888, 1083, 105
Tako, 1844, 1890, 1180, 525
Ebi, 560, 7350, 1420, 525
KurumaEbi, 222, 8925, 6508, 2625
Koiwashi, 1316, 1155, 617, 328
ObaIwashi, 2716, 499, 267, 175
Mentai, 678, 1155, 859, 394
Hamachi, 4772, 998, 632, 105
Isaki, 268, 2625, 1465, 840

このデータは多次元データになっている。
このデータにおいては、対象が魚の名前で、属性がそれぞれ、
卸売数量（Kg）：Oroshi_Kg、高値：Taka、中値：Naka、安値：Yasu、
となっている。つまり、四次元のデータとして表現されていることになる。

ところで、データのことを事前に知っておくことは重要である。

このデータの出典は、広島市中央卸売市場の水産物市況（2012/08/18）から
欠損データを抜いたもの。すでに、ある程度要約されているデータである。

高値、中値、安値、というのは魚の値段のおおよその基準。
魚には大小があって、大小によって値段も、可能な調理の幅も異なる。
高値は、大きくて色々な用途に使えるもの。そして、高い。
中値は、一般的なサイズで、多くの人が目にして、買って、家庭で調理するもの。
安値は、小ぶりなもので、利用の用途が限られるようなもの。
このデータでは、（おそらく）それぞれの値の平均。詳細は不明。

この話に関しては、以下のページに解りやすい説明があった。
http://www.shonai-nippo.co.jp/square/feature/food/sf48.html

少し話が逸れたが、まずは、卸売数量の平均を求めてみる。
Rでは、mean()関数で平均を計算できる。非常に、簡単。

mean(X$Oroshi_Kg)

このコマンドでは、Xという変数の属性名を「$」記号を挟んで指定している。
あるいは、次ようにして求めることもできる。

mean(X[,1])

この表記は、Xの「1列目」という考え方で指定している。
ちなみに、大カッコ（「[ , ]」）のカンマの右側は行数、左側は列数を表す。
以下は、実行結果。

> mean(X$Oroshi_Kg)
[1] 1355.276
> mean(X[,1])
[1] 1355.276

同じことを行なっているので、計算結果は同じ。これは当然。
同様に、他の変数の平均も計算してみる。

mean(X$Taka)
mean(X$Naka)
mean(X$Yasu)

以下は、上記コマンドの実行結果。

> mean(X$Taka)
[1] 3035.103
> mean(X$Naka)
[1] 1453.448
> mean(X$Yasu)
[1] 458.9655

ついでなので、以下のコマンドも実行してみる。
なお、大文字と小文字に注意すること。
Rのコマンドは大文字と小文字を区別する。

colMeans(X)

どのような結果が出てくるか？以下は実行結果。
4つの変数の平均を同時に計算できる。

> colMeans(X)
Oroshi_Kg       Taka       Naka       Yasu 
1355.2759 3035.1034 1453.4483  458.9655 

このコマンドは「列平均（Column Means）」を計算する関数であり、
colMeans()関数と呼ぶ。複雑な分析を行うようになると頻繁に登場する。
今の段階で覚えておくと良い。

今回のデータでは、全く意味は無いが、「行平均（Row Means）」もある。
これは、rowMeans()関数で計算できる。
あくまで練習。これもやってみる。実行結果は以下の通り。

> rowMeans(X)
      Tai     Chinu    Sawara      Yazu    Suzuki      Akou     Kochi     Okoze 
  1439.00    368.50    971.50   1542.00    966.00   2362.50    803.00   2030.75 
   Ainame      Hage Konoshiro    Sayori     Anago    Mebaru   Tachiuo      Hamo 
  1028.75    922.00    831.50   2609.50   1653.25   1631.25   1139.00   1551.25 
    Karei    Hirame       Aji      Saba       Ika      Tako       Ebi KurumaEbi 
  2511.00   1086.50   2324.50   2441.50   1622.50   1359.75   2463.75   4570.00 
 Koiwashi ObaIwashi    Mentai   Hamachi     Isaki 
   854.00    914.25    771.50   1626.75   1299.50 

となる。各対象（魚）について、卸売数量、高値、中値、安値、の平均が出ている。
全く意味は無い。このような無意味なことはやってはいけない。

そもそも、卸売数量と、他の3つの属性は単位が異なっていて比較できないし、
高値、中値、安値というのは、各対象の何かの代表的な値であって、
しかも、中値がある意味でデータの「真ん中」を示している。
あくまで、コマンドの実行例として割りきる。

さて、今回は、「平均」というデータの真ん中について考えた。
ふむ。納得できた。だが、ここで一つの疑問が湧いてくる。
すなわち、平均以外の方法はあり得るのか？という疑問である。

なるほど。確かに、平均という考え方は的を得ているように思えるが、
何か、重大な問題を見落としているようにも思える。何が問題なのか？
この疑問については、次回に考えてみることにしよう。

文系のための「擬逆行列」

さて、前回の話の最後に、「逆行列」というをやった。
「逆行列」というのは、スカラーの「逆数」に対応し、
元の行列に、逆行列を掛けると、「単位行列」が出てくるのであった。

覚えていない人は、もう一度、逆行列の投稿を参照すること。
この投稿の最後に、「逆行列」が「正方行列」でないといけない、と言った。

ところが、この「正方行列」という制限は、色々と不便なのである。
しかし、方法が無いかと言われると、そうでもない。
近似的に「逆行列」を求める方法がある。
っと、いうことを今日は書いてみる。

では、早速、「R」を使って、この問題に挑戦してみる。
まずは、4行×3列の行列を準備する。つまり、n=4、p=3 の行列である。
「R」を立ち上げたら、次のように入力する。

A <- matrix(c(50,0.5,50,5,7,0.2,0,6,75,3,20,50), nrow=4, ncol=3)

今回は、語呂でもなく、特に意味の無いデータであるが、
とにかく、この行列をつかって、演習を行う。

何のエラーも出なければ成功。これで、Aという変数に行列が代入された。
どのようなデータであるのか、変数名を入力してエンターキーを押すと、
変数に格納されているデータが表示される。この例では、変数は「A」である。

> A
     [,1] [,2] [,3]
[1,] 50.0  7.0   75
[2,]  0.5  0.2    3
[3,] 50.0  0.0   20
[4,]  5.0  6.0   50

ふむ。一応、逆行列が計算できるか、確認しみる。
逆行列のコマンドは何であったか？

> solve(A)
Error in solve.default(A) : only square matrices can be inverted

やはり、エラーが出る。「正方行列」でないと計算できない。

困ったことになった。どうしようも無いのか？

仕方が無いので、少し、物の見方というのを変えてみる。

何とか、元の行列を正方行列に変換できないのか？

つまり、という見方である。そのようにすれば、計算できるかもしれない。

ふむ。このような方法は使えるだろうか？

元の行列を転置し、その行列に元の行列を掛けてみる。

説明が解りにくい。これは、文系であっても、数式で示した方が理解しやすい。

そういった事はよくある。完全に数式を避けるべきではない。まぁ、それは良い。

要するに、 ということである。

この式を「R」で計算するには、以下のようにする。

t(A) %*% A

以下は計算結果である。

> t(A) %*% A
        [,1]   [,2]   [,3]
[1,] 5025.25 380.10 5001.5
[2,]  380.10  85.04  825.6
[3,] 5001.50 825.60 8534.0

なるほど、「正方行列」の形になっている。
これを使って「逆行列」を「擬似的」に計算するには、
以下のように計算する。



したがって、「R」では以下のようになる。

> solve(t(A)%*%A)%*%t(A)
            [,1]       [,2]        [,3]        [,4]
[1,]  0.07946341 -0.2671138 -0.04841192 -0.08380351
[2,]  1.54547248 -5.3787729 -1.34597554 -1.45709213
[3,] -0.18729532  0.6772539  0.16092918  0.19593608

ふむ、なるほど。よく見れば、上記の式と対応関係にあることが分かる。

これが「擬似的」に求めた「逆行列」である...本当か？疑いたくなる。
どうすれば、確認ができるのだろうか？確か、何かを掛けると、何かになった。
逆行列が逆数に対応するのだったことを思い出す。

> solve(t(A)%*%A)%*%t(A) %*% A
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -6.095124e-14 -7.318590e-13
[2,] -7.414513e-12  1.000000e+00 -8.988366e-12
[3,]  2.988720e-13 -3.286260e-14  1.000000e+00

つまり、元の行列を掛けると、「単位行列」が出てくるのであった。

なるほど、確かに、対角成分が全て「1」となっていて、

それ以外の成分が全て「0」になっている。逆行列になっている。

このように、擬似的に求めた逆行列のことを「擬逆行列」と呼ぶ。

これは、「正方行列」でなくても使うことができる。これは便利である。

ちなみに、今回の例では、n > p の行列であったが、
p > n の場合には、式が次のようになるので注意が必要。
行列の掛け算は、掛ける方と掛けられる方の順番が重要なのであった。

n > p の場合：
n ＜ p の場合：

実を言うと、「R」では、「擬逆行列」を計算する方法がいくつかある。
これを利用するには、「MASS」というパッケージを読み込む方法がある。
とにかく、次のようにする。

library(MASS)

これで、MASSというライブラリを読み込むことができた。

「R」というソフトウェアは、このようにライブラリと呼ばれるものを読み込み、

様々な機能を追加することができるのである。

このMASSというライブラリを読み込むと擬逆行列を計算する関数が使える。

その関数とは、ginv()関数である。とにかく、使ってみる。

ginv(A)

実行すると、次のようになる。

> ginv(A)

            [,1]       [,2]        [,3]        [,4]

[1,]  0.07946341 -0.2671138 -0.04841192 -0.08380351

[2,]  1.54547248 -5.3787729 -1.34597554 -1.45709213

[3,] -0.18729532  0.6772539  0.16092918  0.19593608

先に、手計算した結果と見比べてみると、なるほど、同じである。

回りくどいことせずとも、逆行列が計算できる。便利である。これを使おう。

この計算が、内部的にされているのだろう...

と思ってはいけない。実は、少々、手の混んだことをやっている。

実は、「擬逆行列」というのは、「特異値分解」からも導くことができ、

ginv()関数では、「特異値分解」を使って計算しているのである。

まだ、「特異値分解」の説明を見ていない人は、

とにかく、先にそっちを確認すること。

さて、特異値分解を行うと、ある行列 X は次のように分解されるのであった。

Uは、「左特異ベクトル」と呼ばれる行列であり、
Dは、「特異値行列」と呼ばれる対角行列。そして、
V'は、「右特異ベクトル」と呼ばれる行列、である。

あまり、深く考えてはいけない。スカラーの「因数分解」のようなものである。
そういう分解をしているのである。つまり、無秩序な任意の行列を、
ある「規則的な行列同士」の「掛け算」に分解しているのである。
複雑な事象を、本質を変えずに、見通しの良い状態にしている、と言っても良い。
ということがニュアンスとして理解できていたら、ここでは良いだろう。

とにかく、この分解された要素の組み合わせ次第で、様々な、難問が解ける。
まさに、数学の魔法が詰まった分解なのである。
擬逆行列の計算というものは、特異値分解の魔法を応用した一例にすぎないのだが、
分解された要素を次のように組み替えれば、擬逆行列となる。



では、「R」で行うとどのようになるのか、やってみる。
まず、変則的な命名方法であるが、A.svd という変数をつくり、
この変数に特異値分解の結果を格納する。

A.svd <- svd(A)

特異値分解は、svd()関数を使えば簡単に計算できる。
また、計算結果がどのように格納されているかは、次の関数で確認できる。

str(A.svd)

つまり、構造（structure）を確認する関数であり、str() というのは、
英語の最初の三文字を取っている。プログラミングでは、このような省略を用いる。
長い変数名は、バグの遠因であるし、また、コーディングも大変である。
実は、変数の命名法にもいくつか流儀というのがあるのだが、この話はいずれ。

すぐに、話が逸れてしまうので困る。私の性格の問題ではあるが。
できれば、愛嬌とか、ユーモアで片付けてもらいたい。まぁ、良いのだが。

つまり、str()関数 の実行結果は次のようになる。


> str(A.svd)
List of 3
 $ d: num [1:3] 110.209 38.706 0.167
 $ u: num [1:4, 1:3] -0.82 -0.0249 -0.409 -0.3995 -0.0647 ...
 $ v: num [1:3, 1:3] -0.5758 -0.0739 -0.8142 0.8161 -0.1118 ...

そして、A.svd という変数に格納されている d, u, v の3つの要素は以下のように取り出す。

> # 特異値行列の抽出

> A.svd$d

[1] 110.2092350  38.7064165   0.1669015

> # 左特異ベクトルの抽出

> A.svd$u

            [,1]        [,2]       [,3]

[1,] -0.82003552 -0.06468901  0.2601641

[2,] -0.02491042 -0.03398276 -0.9059111

[3,] -0.40900495  0.76122378 -0.2263649

[4,] -0.39954495 -0.64435926 -0.2457616

> # 右特異ベクトルの抽出

> A.svd$v

           [,1]       [,2]        [,3]

[1,] -0.5758338  0.8160908  0.04910397

[2,] -0.0738822 -0.1117586  0.99098508

[3,] -0.8142216 -0.5670148 -0.12464897

ふむ。では、まずは、以下のようにしてみよう。

D <- diag(A.svd$d)

V <- A.svd$v

U <- A.svd$u

こうして、3つの変数、D, V, U ができた。

なお、Dに関しては、ベクトルとして出力されているため、

対角行列を作るための diag()関数 を使って、行列に変換している。

さて、ここまで遠回りが続いているが、

公式に基づいて、「擬逆行列」を計算すると以下のようになる。

> V%*%solve(D)%*%t(U)

            [,1]       [,2]        [,3]        [,4]

[1,]  0.07946341 -0.2671138 -0.04841192 -0.08380351

[2,]  1.54547248 -5.3787729 -1.34597554 -1.45709213

[3,] -0.18729532  0.6772539  0.16092918  0.19593608

見事に、擬逆行列が計算されているのが分かる。

特異値分解は、後に主成分分析や対応分析にも登場するので、

今の段階で、どういった特徴を持つ行列であるか、理解しておきたい。

文系のための「逆行列」（1）

スカラーにおける「割り算」は、行列の何に対応するか？
間違っても、スカラーの時のような割り算をしてはいけないし、
分数の表現をするのもNG。ダメダメ。

ふむ、そもそも、「割り算」とはどのようなものだったか？
まずは、そこから整理をしてみよう。

例えば、ケーキを4等分することを考えよう。
1÷4= 0.25 あるいは 1/4 となる。確かにそうなる。

ここで、考え方を少し捻ってみる。いや、逆の考え方をしてみる。
つまり、四等分されたケーキは、いくつ集まって一つになるか？。
これも簡単な問題である。言うまでもなく、「4」である。
「1/4」が「4」個集まって「1」つのケーキになる。

さて、あまりにも当たり前すぎることであるが、
ある数があって、その数を「単位1個分」としたとき、
その数を分母とする分数のことを「逆数」と呼ぶ。

したがって、以下のように言うことができる。

「ある数」 × 「ある数の逆数」 = 「単位数」

ここで、単位数は、例外無く「1」である。
ある数を「単位1個分」とするので当然である。
当たり前すぎるので、文系頭脳では益々混乱の原因になるのだが...。

さて、話を行列に戻そう。つまり、スカラーの逆数に相当するのは何か？
これこそが、「逆行列」と呼ばれるものである。
つまり、ある行列があって、その行列を一つの単位とするとき、
スカラーの単位数に相当するのが「単位行列」であり、
「逆数」に相当するのが、「逆行列」なのである。

かなり、乱暴な説明なので、専門家からは苦言が出てくるかもしれない。
でも、まぁ、数学の知識が皆無であるような、
そういう、自称「文系」向けの話なので、目をつむってもらいましょう。

さて、ここまでで、全体の様子は理解してもらえたとしよう。
ここからが、少々、厄介な話となってくる。

行列の場合、「単位1個分」つまり「単位行列」とは何か、という問題がある。

単位行列は、元の行列に掛けても、
元の行列が出てくるだけ、
という原理を満たさないといけない。

そのようになるには、対角成分が「1」でそれ以外が「0」になる行列しか無い。
ためしに、適当な例をやってみる。



全部を計算するのは大変なので、一行目のだけ手計算。
行列の「かけ算」の仕方を実践すると、

掛けられる方の行列は「横方向」に1, 9, 4 となり、
掛ける方の行列は「縦方向」に1, 0, 0 となる。
それで計算すると、上のようになる。これが解の一行目。
なるほど、確かに元の行列と同じになるようだ。

最後にまとめに入ろう。要するに、「逆行列」とは何か？
我々が通常「割り算」と読んでいるものは、観点を変えると、
スカラーの世界においては「逆数」を掛けるということであり、
行列の世界においては「逆行列」を掛けるということなのである。
つまりは、「ある行列」に「別の行列」を掛けて「単位行列」となるとき、
その「別の行列」のことを「逆行列」と呼び、慣例的には、
行列の変数の右上に「-1」を付けて、 と表す。

一般的に、「逆行列」というのは、「正方行列」にしか適応できず、
縦と横の長さが異なる場合には、「疑逆行列」というのが必要になる。

文系のための「特異値分解」（２）

ふむふむ。何となく、「特異値分解」というものが、
複雑な状況を整理していることが解ってもらった...と信じたい。
要するに、難しく考えてはいけないのである。肩の力を少し抜いてみよう。

っで、肩の力を抜いたところで、
今度は、肩に力が入ってしまいそうな話をする。
ここでは、ちょっと我慢というのが必要である。

つまり、「特異値分解」がどういった分解であるか？ということである。
簡単に言うと、あるn行 × p列 から成る行列があって、
その行列を三つの行列に分解することである。すなわち、

• 一つ目の行列が、「左特異ベクトル」と呼ばれる行列（u）
• 二つ目の行列が、「特異値行列」と呼ばれる行列（d）
• 三つ目の行列が、「右特異ベクトル」と呼ばれる行列（v）

この三つのうち、「u」と「v」が「直交行列」と呼ばれる行列であり、
「d」は「対角行列」と呼ばれる行列で、
その対角成分に「特異値」というのが並んで入っている。

「特異値行列」とは、全くもって、不可解な用語であるが、
言葉の意味から入ると、混乱の原因となるので、
気持ち悪くても、以下のような性質があるということで納得したい。

まず、「特異値行列」である「d」というのは、
ある次元を、別の次元に変換するためのある種の変換行列であって、
必ず、正の値を取るようになっている。

二つの「特異ベクトル」行列は、「直交行列」であることが重要。
というのは、「直交行列」は定義上、その転置行列が「逆行列」に等しく、
つまり、元の「直交行列」にその転置行列を掛けると「単位行列」が出てくる。

この性質を用いることで、様々な複雑な計算が魔法のように解ける。
「対象の数（n）」が「変数の数（p）」よりも多い状況を想定すると、
u の行列は元の行列と同じサイズの行列で、d と v の行列はp× pの正方行列になる。

こういう分解をするのが「特異値分解」と呼ばれる分解である。

では、なぜ、このような不思議な変換を用いるかと言うと、
この分解を行うことで「擬逆行列」を求めるのが非常に楽になるほか、
「主成分分析」の計算においても利用される。