金子邦彦研究室プログラミング主成分分析,次元削減R システムを用いた不偏分散共分散行列

R システムを用いた不偏分散共分散行列

このページでは,R システムを用いた不偏分散共分散行列の求め方と図解で説明するとともに, 不偏分散共分散行列の意味を説明する. (ここでは説明しませんが octave の cov() 関数を使った場合も,R と同様の結果が出てきます).

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

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

[image]
[image]


前準備

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

R システムの CRAN の URL: https://cran.r-project.org/

CSVファイルを読み込み,データフレームに格納

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

    Book1.csv をダウンロード (参考: 「外国為替データ(時系列データ)の情報源の紹介」の Web ページ

    以下の説明では、

    として説明を続ける.

    自前の CSV ファイルを使うときの注意: read.table() 関数を使うので, 属性名は英語になっていること.属性名は,CSV ファイルの第一行目に書いていること.

  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);
    

    R の read.table のオプション

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

    次のコマンドを実行.

    edit(X);
    

    次のコマンドを実行.

    str(X)
    

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

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

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

不偏分散共分散行列が結果として得られる. 結果が3行3列の配列になっている

cov( X[,3:5] )

[image]

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

概要

R の cov() 関数を使う方法だと、中身がブラックボックスです. 実は,以下の 2 つの R の式が,同じ結果になります. ここでは、確かに同じ結果になることを確認する手順を,図解などで説明する.

[image]

手順

上では「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);
    

    [image]

    ◆ Linux での動作手順例

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

    [image]

    read.table() のオプション

  2. オブジェクト X の確認

    次のコマンドを実行.

    edit(X);
    

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

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

    colMeans( X[,3:5] )
    

    [image]
  4. 行列 X の各要素から,X の列の平均値を引く

    つまり、次のことを行いたいのです.

    ※ ここで行っていることは, 変数 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 ) 
    

    実行結果の例は次の通りです.

    [image]

    今度は,要素 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 ) 
    

    実行結果の例は次の通りです.

    [image]

    行列の行数と,ベクトルの長さが等しいときは, 行列の各要素から,ベクトルの値を引く(行列の行番号と、ベクトルの要素番号が等しい)ことが見て取れます.

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

    ※ データを転置したいので関数 t() を使う.

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

    [image]

    上のように書く代わりに、次のように書いても同じ意味ですが、リサイクル規則を使わないと、書くのが面倒です ※ 「V[,3] - colMeans(V[,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 は、610行3列の行列になっていることに注意.