Rを用いた不偏分散共分散行列

第1章　不偏分散共分散行列と不偏相関係数行列

3 つの変数を n 回観測して得られる n 個の標本から構成される長さ n のベクトル に対する不偏分散共分散行列と不偏相関係数行列の定義は次の通り．

• 分散には，標本分散と不偏分散の2種類がある．同様に， 共分散にも，標本分散の拡張と不偏分散の拡張の2種類がある． R や octave では，不偏分散の方が使われる（標本分散ではない）．
• ここでは3 つの変数としているが，一般の m 個の変数が与えられた場合も同様の定義になる．

第2章　前準備

R システムのインストール

• Windows での R システムのインストール: 別ページ »で説明
• Ubuntu での R システムのインストール: 別ページ »で説明

第3章　CSVファイルを読み込み，データフレームに格納

1. (前準備) 使用する CSV ファイルの作成

   Book1.csv をダウンロードする （参考： 「外国為替データ（時系列データ）の情報源の紹介」の Web ページ）．

   以下の説明では，

   • Windows の場合：　データファイル名： C:\R\Book1.csv
   • Linuxの場合：　データファイル名： /tmp/Book1.csv

   として説明を続ける．

   * 自前の CSV ファイルを使うときの注意： read.table() 関数を使うので， 属性名は英語にする．属性名は，CSV ファイルの1行目に書く．
2. 使用する CSV ファイルの確認

   属性名が CSV ファイルの1行目に書かれていることを確認する．

   Book1.csv
3. R の起動
4. read.table 関数を用いて，CSV ファイルを R のデータフレームに読み込む

   次のコマンドを実行する．

   ◆ Windows での動作手順例

   X <- read.table("C:/R/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
   

   ◆ Linux での動作手順例

   X <- read.table("/tmp/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
   

   read.table 関数のオプション

   • X <- ・・・ 変数 X に読み込むという意味
   • C:/R/Book1.csv, "/tmp/Book1.csv" ・・・ 読み込む CSV ファイル名．Windows では区切りに「/」を使うことに注意．
   • header="TRUE" または header="FALSE" ・・・ 列ラベルが設定されているか
   • sep="," や sep="\t" や sep=" " など ・・・ 列を区切る記号（CSV ファイルのときは「sep=","」）
   • na.strings="NA" ・・・ 欠損値には "NA" を使うという意味
   • dec="." 　・・・　ファイルで使われている小数点記号（既定値はピリオド）
   • strip.white=TRUE ・・・ 個々のデータの先頭や末尾にある空白文字を取り除いて読み込む
   • skip=＜行数＞ ・・・ 読み飛ばし行数
   • nrow=＜行数＞ ・・・ 読み込み行数
5. オブジェクト X の確認

   次のコマンドを実行する．

   edit(X);
   

   次のコマンドを実行する．

   str(X)
   

第4章　cov() 関数を使って不偏分散共分散行列を求める手順

cov() 関数を使って不偏分散共分散行列を求める

次のコマンドを実行する． データはUSDとEURとAUDの3列（それぞれ3列目と4列目と5列目）なので， 「3:5」と書く．これは「c(3,4,5)」と同じ意味である．

不偏分散共分散行列が結果として得られる． 結果は3行3列の行列になる．

cov( X[,3:5] )

第5章　より基本的な関数を使って不偏分散共分散行列を求める手順

概要

R の cov() 関数を使う方法では，計算の中身が見えない． 以下の2つの R の式は同じ結果になる． ここでは，同じ結果になることを図解で確認する．

• cov() 関数を使って不偏分散共分散行列を求める R の式

  cov( X[,3:5] )
  

• colMeans() 関数，t() 関数などを使って不偏分散共分散行列を求める R の式

  crossprod( D ) / ( nrow(D) - 1 )
  * この式は「t( D ) %*% D / ( nrow(D) - 1 )」と同じ意味
  

  但し，D は次の手順で作る． D は，行列 X の各要素から各列の平均値を引き，転置した行列である． 「転置」にしているのは，R のリサイクル規則を使って簡単に作れるためである．

  D <- t( t( X[,3:5] )  - colMeans( X[,3:5] ) )
  

  この Web ページの先頭で説明したように，不偏分散共分散行列は，不偏分散を拡張した形の分散共分散行列になっている．

手順

上では「crossprod( D ) / ( nrow(D) - 1 )」， 「D <- t( t( X[,3:5] ) - colMeans( X[,3:5] ) )」のように書いた． 式が複雑なので，一歩ずつ確かめながら R の理解を深める．

1. read.table 関数を用いて，CSV ファイルを R のデータフレームに読み込む

   次のコマンドを実行する．

   ◆ Windows での動作手順例

   X <- read.table("/tmp/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
   

   

   ◆ Linux での動作手順例

   X <- read.table("/tmp/Book1.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE);
   

   

   read.table 関数のオプション

   • X <- ・・・ 変数 X に読み込むという意味
   • C:/R/Book1.csv, "/tmp/Book1.csv" ・・・ 読み込む CSV ファイル名．Windows では区切りに「/」を使うことに注意．
   • header="TRUE" または header="FALSE" ・・・ 列ラベルが設定されているか
   • sep="," や sep="\t" や sep=" " など ・・・ 列を区切る記号（CSV ファイルのときは「sep=","」）
   • na.strings="NA" ・・・ 欠損値には "NA" を使うという意味
   • dec="." 　・・・　ファイルで使われている小数点記号（既定値はピリオド）
   • strip.white=TRUE ・・・ 個々のデータの先頭や末尾にある空白文字を取り除いて読み込む
   • skip=＜行数＞ ・・・ 読み飛ばし行数
   • nrow=＜行数＞ ・・・ 読み込み行数
2. オブジェクト X の確認

   次のコマンドを実行する．

   edit(X);
   

   
3. X の各列の平均を求める

   データはUSDとEURとAUDの3列（それぞれ3列目と4列目と5列目）ある．この3列について， 各列ごとに平均を求める． 「3:5」は「c(3,4,5)」と同じ意味である．平均を求めた結果はベクトルになる．

   colMeans( X[,3:5] )
   

   
4. 行列 X の各要素から，X の列の平均値を引く

   つまり，次のことを行いたい．

   • X の3列目の要素すべてから，colMeans( X[,3:5] ) の第1要素を引く
   • X の4列目の要素すべてから，colMeans( X[,3:5] ) の第2要素を引く
   • X の5列目の要素すべてから，colMeans( X[,3:5] ) の第3要素を引く

   * ここで行っているのは， 変数 X 内の3列目と4列目と5列目がなす「3次元のベクトル集合」の平均が原点になるように平行移動させる操作とみなせる．

   行列からベクトルを引くとき，リサイクル規則が使われる． ここで，リサイクル規則を例を使って説明する． 次のような，要素 1, 2, 3, 4, 5, 6 が入った3行2列の行列があるとする．

   1 4
   2 5
   3 6
   

   この行列に対してベクトル c(2, 5) の引き算を行うと，リサイクル規則により， ベクトル c(2, 5) が次のような行列に引き伸ばされる．

   2 5
   5 2
   2 5
   

   matrix( c( 1, 2, 3, 4, 5, 6 ), nrow = 3, ncol = 2 )
   matrix( c( 1, 2, 3, 4, 5, 6 ), nrow = 3, ncol = 2 ) - c( 2, 5 ) 
   

   実行結果の例は次の通りである．

   

   今度は，要素 1, 4, 2, 5, 3, 6 が入った2行3列の行列があるとする．

   1 2 3
   4 5 6
   

   この行列に対してベクトル c(2, 5) の引き算を行うと，リサイクル規則により， ベクトル c(2, 5) が次のような行列に引き伸ばされる．

   2 2 2
   5 5 5
   

   matrix( c( 1, 4, 2, 5, 3, 6 ), nrow = 2, ncol = 3 )
   matrix( c( 1, 4, 2, 5, 3, 6 ), nrow = 2, ncol = 3 ) - c( 2, 5 ) 
   

   実行結果の例は次の通りである．

   

   行列の行数とベクトルの長さが等しいときは， 行列の各要素から，行番号に対応するベクトルの値を引くことが分かる．

   「行列の各要素から，列の平均値を引く」ことをしたいので， 行列 X[,3:5] を転置して2834行3列の行列を作り，列の平均値のベクトル（長さ3）で引き算する． リサイクル規則が使えるので，次のような簡単な式になる（結果は2834行3列の行列である）．

   * データを転置するには関数 t() を使う．

   D <- t( t( X[,3:5] )  - colMeans( X[,3:5] ) )
   edit(D); 
   

   

   上のように書く代わりに，次のように書いても同じ意味であるが，リサイクル規則を使わないと記述が長くなる． * 「X[,3] - colMeans(X[,3])」は，3列目の全要素から3列目の平均値を引くという意味．

   hoge <- matrix(
             c( ( X[,3] - colMeans(X[,3]) ), 
                ( X[,4] - colMeans(X[,4]) ),
                ( X[,5] - colMeans(X[,5]) ) ), nrow = 2834, ncol = 3 ); 
   

5. 不偏分散共分散行列を求める

   crossprod() 関数はクロス積（この場合 t(D) %*% D）を求める関数である．D は2834行3列の行列であることに注意する．

   crossprod( D ) / ( nrow(D) - 1 )