3年生 数理実習Ａ (統計手法) 資料
R
による多変量データ解析
—
グラフィカルモデリング
—
数理3研助教 松井 千尋1 61号講義室, 2016/4/7(木), 4/14(木), 4/21(木)目 次
1 実習の目的 1 2 実習の手順 1 2.1 Rのインストール,基本的な使い方 (本年度より省略) . . . . 2 2.2 パッケージggm のインストール. . . . 3 2.3 サンプルデータの読み込みと,基本統計量による検討 . . . . 4 2.4 グラフの用語 . . . . 6 2.5 グラフィカルモデリングの準備1 (偏相関行列) . . . . 7 2.6 データの収集 . . . . 8 2.7 グラフィカルモデリングの準備2 (無向独立グラフ) . . . . 8 2.8 共分散選択アルゴリズム . . . . 10 2.9 共同実験者のデータの解析 . . . . 14 3 レポート提出要領 141
実習の目的
本実習では,実データの解析を通して,統計学の手法のいくつかを学んでもらう. 特に,３つ以上 の変数の「絡み」を視覚的にとらえる手法として有用な, グラフィカルモデリングについて学ぶ. また,フリーの統計解析ソフト「Ｒ」をインストールし,その環境に慣れることも目的とする. 最 終目標は，集めた実データを学生同士で交換し，お互いのデータを解析することである．2
実習の手順
実習の手順を説明する. 最終的な目標は,自分で収集した統計データから,図2 (p. 10)のような グラフを作ることである. 1 まつい ちひろ, 6号館356号室, matsui@mist.i.u-tokyo.ac.jp; TA荒木 健司. http://www.sr3.t.u-tokyo.ac.jp/matsui/2016s/2016s jikken3.html2.1
R
のインストール, 基本的な使い方 (本年度より省略)
まず, http://cran.md.tsukuba.ac.jp/にアクセスし,「Download R for Windows」→「base」 をクリックする. Rの最新版 (2014-4-17現在) である R-3.1.0-win.exeをダウンロードし, 実 行する2. インストールが終わると, デスクトップにアイコンが作成される. ダブルクリックして,以下の ような画面が現れることを確認する. とりあえず,次のような簡単な計算ができることを確認してみる. Rの基本操作 > 3*5 # かけ算 [1] 15 > c(19,76)+c(11,13) # ベクトルの計算 [1] 30 89 > x <- c(3,1,4,1,5,9) # 代入 > 1:4 # 数列の生成 [1] 1 2 3 4 > 1:4+1 # (1:4)+1 と同じ [1] 2 3 4 5 > f <- function(x,y){ sqrt(x^2+y^2) } # 関数の定義 > f(3,4) [1] 5 Rのコマンドをテキストファイルに羅列してまとめて実行したいときは,以下のようにする. 1. 「ファイル」→「ディレクトリの変更...」として,ファイルのあるディレクトリに移動. 2. 「ファイル」→「Rコードのソースを読み込み...」として,ファイルを選択. プログラムの特定の箇所で実行を止めたいときは,関数 browser() を挿入しておけばよい. その他, Rの使い方全般については,参考文献あるいは下記URLを参照すること. http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html 2既に自分のPCにRがインストールされている場合はこれらの手順を省く．
問 1. Rに付属の以下の関数や演算について，Rの参考書・ヘルプ・webサイトを参考にしながら
動作を確認せよ3．
総和，総積，累積和 sum, prod, cumsum 絶対値，対数，平方根 abs, log, sqrt
平均，標準偏差，最小値，中央値，最大値 mean, sd, min, median, max 行列の各列の標準化 scale ベクトル生成（combine）, ベクトル化 c リスト list 行列, テンソルの生成 matrix, array 行列積 %*% 行列の転置 t 行列の結合 rbind, cbind 対角行列，対角要素，単位行列 diag 相関行列 cor 行列式 det 各行 or 各列の関数評価 apply 連立方程式，逆行列 solve データファイルの読み込み，書き出し read.table, write.table 散布図，対散布図 plot, pairs 図の EPS ファイルへの出力 dev.copy2eps 論理値が真であるインデックスを返す which 正規分布の累積分布，乱数発生 pnorm, rnorm 変数の表示，変数の削除 ls, rm デバッグ browser
繰り返し文 for, while, repeat
2.2
パッケージ ggm のインストール
後でグラフィカルモデリングによる解析を行なう際,さらにパッケージggm (Graphical Gaussian Models) のインストールが必要となる. 以下の手順でインストールする： 1. Rのメニューから「パッケージ」→「パッケージのインストール」を選ぶ (例えば「Japan (Tsukuba)」を選ぶ). 2. 「ggm」を選ぶ.3. “package ‘ggm’ successfully unpacked and MD5 sums checked”などの表示が出ればインス
トール終了. これだけでは（なぜか）依存パッケージであるgraph とRBGL がインストールされない．仕方な いので，手動でインストールする．まずはgraphをRのウィンドウで以下のコマンドを入力して インストールする: graph のインストール > source("http://bioconductor.org/biocLite.R") > biocLite("graph") 同様に RBGL もインストールする。 3問に対する答えはレポートには書かなくてよい．以下同様．
RBGL のインストール > source("http://bioconductor.org/biocLite.R") > biocLite("RBGL") インストールが終わってもそのままでは利用することはできない. 下記コマンドでパッケージを読 み込む（何も出力されないが,それでよい）. このコマンドは, Rを起動するたびに1回実行する必要がある. ただし, テキストファイルの最初の行にこのコマンドを書いておき, 前節の方法で読み込めば, 手 間は省ける. パッケージ読み込み > library(ggm)
2.3
サンプルデータの読み込みと, 基本統計量による検討
ある集団を複数の調査項目に着目して観察し, それをまとめたものを統計データとよぶ. また, 調査項目のことを変数といい,数値で表される変数を量的変数,数値で表されない変数を質的変数 という. 本実習では量的変数だけからなるデータを扱う. 例として, ggmに付属のサンプルデータ marksを使い,解析の流れを見る. データは次のように して読み込む. データの読み込み > data(marks) データの内容は次のようにして見ることができる： データの表示 > marksmechanics vectors algebra analysis statistics
1 77 82 67 67 81 2 63 78 80 70 81 3 75 73 71 66 81 4 55 72 63 70 68 ... (中略) ... 87 5 26 15 20 20 88 0 40 21 9 14 このデータは88人の学生の5教科の成績をまとめたものであり,変数は5個,データ数は88個 と考える. 以後,このデータのことを成績データと呼ぶ. 数式を扱う際は,データを行列と見なし て(Xt,i)1≤t≤n,1≤i≤p (n = 88, p = 5)とおく. 例えばX2,3= 80である. 基本統計量の算出は，関数summary で実行できる．各変数の要約として最小値 (min),第1四
分位点 (1st Qu.),中央値 =第2四分位点 (median),平均値 (mean), 第3四分位点(3rd Qu.), 最 大値 (max) がそれぞれ表示される. ただし,平均値は ¯ Xi= 1 n n ∑ t=1 Xt,i
で定義される量であり,四分位点とはデータを小さい順に並べたときに頭から25%, 50%, 75% の
位置にあるデータのことである.
成績データの要約
> summary(marks)
mechanics vectors algebra analysis statistics Min. : 0.00 Min. : 9.00 Min. :15.00 Min. : 9.00 Min. : 9.00 1st Qu.:30.00 1st Qu.:42.00 1st Qu.:45.00 1st Qu.:35.75 1st Qu.:31.00 Median :41.50 Median :51.00 Median :50.00 Median :49.00 Median :40.00 Mean :38.97 Mean :50.59 Mean :50.60 Mean :46.68 Mean :42.31 3rd Qu.:49.25 3rd Qu.:60.00 3rd Qu.:57.25 3rd Qu.:57.00 3rd Qu.:51.50 Max. :77.00 Max. :82.00 Max. :80.00 Max. :70.00 Max. :81.00
次に相関行列を求めてみる4. 相関行列 (ri,j) は次式で定義される： ri,j = si,j √s i,isj,j , si,j = 1 n n ∑ t=1 (Xt,i− ¯Xi)(Xt,j − ¯Xj) 相関行列の計算 > cor(marks)
mechanics vectors algebra analysis statistics mechanics 1.0000000 0.5526975 0.5462281 0.4096365 0.3894430 vectors 0.5526975 1.0000000 0.6096447 0.4850813 0.4364487 algebra 0.5462281 0.6096447 1.0000000 0.7108059 0.6647357 analysis 0.4096365 0.4850813 0.7108059 1.0000000 0.6071743 statistics 0.3894430 0.4364487 0.6647357 0.6071743 1.0000000 > round(cor(marks),3) # 結果を小数点以下第3桁までに丸めて表示したいとき
mechanics vectors algebra analysis statistics mechanics 1.000 0.553 0.546 0.410 0.389 vectors 0.553 1.000 0.610 0.485 0.436 algebra 0.546 0.610 1.000 0.711 0.665 analysis 0.410 0.485 0.711 1.000 0.607 statistics 0.389 0.436 0.665 0.607 1.000 相関行列は添字iとjに関して対称であることに注意する．これらの量的な要約値に加えて,ヒ ストグラム,散布図といった図的要約を行なうことも重要である. R には,そのようなさまざまな 組み込み関数が用意されている. 例として,各変数の組合せに対する散布図(対散布図) を図1 に 示す. 相関行列や対散布図からサンプルデータ marksに関して分かることとして,全ての教科の間に 正の相関があること,そしてmechanics とstatisticsとの相関は比較的小さいことなどが挙げ られる. このように,相関行列や対散布図は2つの変数間の関連性を見るのには適しているが，３ つ以上の変数が互いにどう絡み合っているかについては,あまり直接的な示唆を与えない. 例えば,
vectorsの成績を固定した下でmechanics とstatisticsの相関はあるか,などという情報は分
かりづらい.
mechanics 20 60 10 30 50 70 0 20 60 20 60 vectors algebra 20 40 60 80 10 30 50 70 analysis 020 60 20 40 60 80 10 40 70 10 40 70 statistics 図1: 対散布図. > pairs(marks) このように, 2つの変数を比較したいときに他の変数を固定するという方法は, 科学的実験では よく用いられる. 例えば,電気抵抗の測定において,温度一定のもとで電圧と電流を測る,といった 具合である. しかし,成績のようなデータでは,他の変数を固定して2変数間の関係を観察すると いうわけにはいかない. 次節以降で扱うグラフィカルモデリングは，この問題を解決してくれる．
2.4
グラフの用語
グラフGとは,頂点集合V と辺集合E⊂ V ×V の組で定義されるオブジェクトである. 本実習で は,辺の向きを考慮しないグラフ(=無向グラフG = (V, E))のみを扱う. 例えばV ={a, b, c, d}, E = {(a, b), (a, c), (b, d)}に対応するグラフをggm で描きたいときは以下のコマンドを実行する.なお, drawGraph の引き数であるadjustをTRUE にすると,マウスを使ってグラフを調整できる.
グラフの調整は右クリックで停止できる.
グラフの描画
> amat <- UG(~ a*b + a*c + b*d) > amat a b c d a 0 1 1 0 b 1 0 0 1 c 1 0 0 0 d 0 1 0 0 > drawGraph(amat,adjust=FALSE) a b c d
2.5
グラフィカルモデリングの準備 1 (偏相関行列)
グラフィカルモデリングとは,一言でいえば,データの従うシンプルなモデルを相関行列の逆行 列に関する制約から選び(共分散選択), 選択された関連構造を,無向グラフの形で表すものである. グラフィカルモデリングの理論の詳細は, webページの付録および参考文献の宮川(1997)などに 譲り,ここでは実データへの適用方法に重点をおいて説明する. いま,相関行列をR = (ri,j) と表し,その逆行列を R−1 = (ri,j) と表記する. このとき, pi,j = −r i,j √ ri,i√rj,j (i̸= j), 1 (i = j) (1) は,変数i と 変数j の偏相関係数と呼ばれる. 偏相関係数を並べた行列 P = (pi,j) が偏相関行列 である. 成績データに関して, 相関行列 (再掲) と偏相関行列を計算したものを表 1 に示す. 両者 を区別するために,ここでは偏相関行列の対角成分は — で表すという慣習にしたがう. また, と もに対称行列であるため, 右上成分は表示していない(レポートではこの記法に従う必要はない). さらに,変数名をS1 = mechanics, S2 = vectors, S3 = algebra, S4 = analysis, S5 = statistics
と略記した. 表1: 成績データの相関行列 (左) と偏相関行列 (右). S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.41 0.485 0.711 1 S5 0.389 0.436 0.665 0.607 1 S1 S2 S3 S4 S5 S1 —— S2 0.328 —— S3 0.229 0.282 —— S4 -0.001 0.078 0.432 —— S5 0.026 0.02 0.357 0.253 —— 問 2. 以下の2つの関数はともに, 相関行列R を入力とし, 偏相関行列P を出力する関数である ことを確かめよ. 偏相関行列を計算する関数I cor2par <- function(R){ X <- solve(R) p <- nrow(R) P <- matrix(0,p,p) dimnames(P) <- dimnames(R) for(i in (1:p)){ for(j in (1:p)){
if(i != j) P[i,j] <- -X[i,j]/sqrt(X[i,i]*X[j,j]) if(i == j) P[i,j] <- 1
}} P }
偏相関行列を計算する関数II cor2par <- function(R){ X <- solve(R) d <- sqrt(diag(X)) P <- -X / (d %*% t(d)) diag(P) <- 1 P }
2.6
データの収集
ここからは，ここまでで学んだ方法を使い，各自で集めてもらったデータを解析してもらう．ど のようなデータでもよいが,変数の個数は 4∼ 8 程度が望ましい. 例えば, 都道府県別に人口, 平 均所得,進学率,就職率を変数とした47× 4のデータなどが考えられる. また 総務省統計局 http://e-stat.go.jp/ StatLib http://lib.stat.cmu.edu/datasets/ などの webサイトも参考になる. 自前のデータの読み込みは以下の手順で行う．まず，集めた以 下のようなデータ:"math" "phys" "chem" "eng" "1" 30 40 50 60 "2" 20 0 10 70 "3" 40 60 90 45 "4" 30 60 20 70 "5" 40 50 20 70 "6" 50 70 30 80 をファイル名mark4.txtとして保存する．Rのメニューから「ファイル」→「ディレクトリの変 更...」として, データファイルのあるディレクトリに移動し,その後次のコマンドで自前データの 読み込みができる: 自前のデータの読み込み > X <- read.table("mark4.txt") 問3. インターネットや年表などを使ってデータを収集し,各基本統計量（要約・相関・偏相関）の 計算や散布図の描画を試してみよ. ここで集めたデータは2.9 節で用いるので保存しておくこと. 数値データに含まれるケタ区切りのカンマは不具合の原因となるので取り除くこと(例えばWord などの置換を用いよ．) （実習１日目はここまでを目標とする）
2.7
グラフィカルモデリングの準備 2 (無向独立グラフ)
変数iと変数 j の偏相関係数 pi,j の意味は,「変数 iと変数 j 以外のすべての変数を固定したときの,変数iと変数j の相関の指標」 である. 正確には次の定理のようにまとめられる. 定理 1 (多変量正規分布における偏相関係数). データが多変量正規分布に従っていると仮定する. このとき, 変数 i と変数j は, その他のすべての変数を与えた下で,相関 pi,j の正規分布に従う. 特に,変数 i と変数j が条件付き独立であるための必要十分条件はpi,j = 0 である. したがって, pi,j = 0ならば,変数 iと変数 jは,残りの変数を固定したときに無相関になる. 成 績データでは,例えば, analysis とmechanics の偏相関係数(p41=−0.001)は値が 0 に非常に近 い. これは,「もし, vectors, algebra, statistics全ての成績が同じ人だけを集めることができれば,
彼らの analysisとmechanics の成績に相関はほとんどない」ことを意味する. グラフィカルモデルとは,偏相関行列のうちゼロに近い要素は積極的にゼロとおくことにより, データの発生メカニズムをより簡潔に表したモデルである5. このとき, • 各変数をグラフの頂点に対応させる. • 偏相関行列が非ゼロの変数対に対応する頂点間に辺を引く. • 偏相関行列がゼロの成分に対応する2頂点間には辺を引かない. として作られるグラフを 無向独立グラフ と呼ぶ. 表1の偏相関行列に対して, 次節で述べるアルゴリズムによって共分散選択を行なったところ, 10個の偏相関係数のうち4個をゼロとおいたモデルが選択された. 選択されたモデルでの相関行 列,偏相関行列の推定値を,表2に示す. ゼロとおいた偏相関係数には下線をつけてある. また,対 応する無向独立グラフを図2 に示す. 表2の相関行列の推定値を表 1の相関行列と比べてみると, 下線をつけている変数対のみ,値がずれていることが確認できる. 表 2: 選択されたモデルでの相関行列の推定値 (左)と偏相関行列の推定値 (右) S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.388 0.433 0.711 1 S5 0.363 0.405 0.665 0.607 1 S1 S2 S3 S4 S5 S1 —— S2 0.331 —— S3 0.235 0.327 —— S4 0 0 0.451 —— S5 0 0 0.364 0.256 —— 次の定理は,グラフィカルモデリングの有効性を示すものである. 証明は宮川(1997)の定理(参 考文献 [2])を参照せよ. 定理 2 (分離定理). 無向独立グラフにおいて, a, b, s を互いに排反な頂点集合とし, aとbを結ぶ 任意の道がsの要素を必ず通るとする(a とbがsで分離されている,という). このとき,変数集 合sを与えたもとで,変数集合aと変数集合b は条件付独立になる. 5 統計学では,確率分布の集合のことをモデルと呼ぶ.本稿では多変量正規分布だけを考えているので,モデルと相関 行列の集合（あるいは偏相関行列の集合）を同一視できる.
図 2: 成績データの無向独立グラフ（偏相関係数は描画ソフトで書き加えた）.
例えば, a = {vectors}, b = {statistics}, s = {algebra}とおいてみると, 図2 より定理 2 の条 件は満たされている. よって, algebra の点数を与えたとき, vectorsと statistics は無相関である ことが分かる(mechanicsとanalysis を与える必要がなくなる).
一方, a ={mechanics,vectors}, b = {analysis,statistics}, s = {algebra} とおいてみると, やは
り定理 2 の条件は満たされている. したがって,例えば mechanics と vectors の関係に直接影響
するのはalgebra のみであることが分かる.
問 4. a ={mechanics}, b = {statistics}, s = {vectors,analysis} とした場合，aと bは sを与え た下で独立となるだろうか？
2.8
共分散選択アルゴリズム
2.7節で結果だけ示した共分散選択 (すなわち,表1 から表2の導出) の方法について説明する. 理論的背景に興味があれば, webページの付録を参照のこと. また,以下の手順を全てまとめて一 つの関数にすることもできるが,対話的に進めていく方が理解が深まると思われるので,あえて冗 長な説明をする. まず,入力・出力を明らかにしておく． 共分散選択アルゴリズムの入力・出力 入力 相関行列R, 標本サイズn (成績データの場合 n = 88). 出力 AICを(なるべく)小さくする無向独立グラフ G, Gに対応する相関行列M = M (G), 偏相関行列P = P (G), および計算途中結果の表示（またはリストにして保存）. ここでいくつか言葉の意味を説明しおく. まず,「無向独立グラフ Gに対応する相関行列」と は,「Gが表す制約の下で,最もデータにフィットするような相関行列」のことであり,具体的には最尤法と呼ばれる方法で求めることになる. また, AIC (Akaike’s Information Criterion) とはモ
デルのあてはまりの良さを表す尺度であり (小さいほど良いモデル),
と定義される. 任意定数は, 制約を置かないモデルのAIC が0 となるように定める (詳細はweb ページの付録を参照). ここでは,パッケージ ggm に用意されている関数 fitConGraph を使って計算する方法を示す. 図2 のグラフに対応する相関行列,およびAIC の導出は以下のプログラミングで行う. 関数 fitConGraph の使い方, AIC の計算法 > options(digits=3) # 表示桁数を3桁までにしたい場合 > X <- marks # データ > n <- nrow(X); p <- ncol(X) # サイズ > R <- cor(X) # 相関 > amat <- matrix(1,p,p)-diag(p); # 全ての頂点対を辺で結んだグラフ, # amat <- UG(~a*b*c*d*e) でも同じ. > dimnames(amat) <- dimnames(R) # 変数名をコピー > amat[4,1] <- amat[1,4] <- 0 # 辺 (4,1) を除去 > amat[4,2] <- amat[2,4] <- 0 # 辺 (4,2) を除去 > amat[5,1] <- amat[1,5] <- 0 # 辺 (5,1) を除去 > amat[5,2] <- amat[2,5] <- 0 # 辺 (5,2) を除去 > amat # 表示
mechanics vectors algebra analysis statistics
mechanics 0 1 1 0 0 vectors 1 0 1 0 0 algebra 1 1 0 1 1 analysis 0 0 1 0 1 statistics 0 0 1 1 0 > f <- fitConGraph(amat,R,n) # 最尤法 > f # 出力結果はリストである $Shat # 推定した相関行列
mechanics vectors algebra analysis statistics mechanics 1.000 0.553 0.546 0.388 0.363 vectors 0.553 1.000 0.610 0.433 0.405 algebra 0.546 0.610 1.000 0.711 0.665 analysis 0.388 0.433 0.711 1.000 0.607 statistics 0.363 0.405 0.665 0.607 1.000 $dev # 逸脱度 [1] 0.9 $df # 制約式の個数 [1] 4 $it # 繰り返し計算における繰り返し数 [1] 2 > f$dev # リストから値を取り出す [1] 0.9
> aic <- f$dev - 2*f$df # AIC > aic
[1] -7.1
共分散選択アルゴリズムは次のように述べられる. 共分散選択アルゴリズム 1. Gを完全グラフ(全ての頂点対を辺で結んだグラフ)とする. M = Rとおく. AICの初 期値を0 とする（完全グラフの AIC）. 2. 相関行列M (R ではない)から偏相関行列 P = (pi,j) を求める(問 2参照).
3. Gの辺(i, j) のうち,偏相関の絶対値 |pi,j|が最小となる(i, j)を選び, Gから取り除く.
4. Gに対応する相関行列M = M (G)とAIC = AIC(G)を求める(fitConGraphを使う．
fitConGraph の第二引数は M でなくR とすること). 5. AICが減少したら手順 2に戻る. AIC が増大したら,その直前のグラフを出力する. 問 5. 上記アルゴリズムの手順3にある「Gの辺(i, j)のうち,偏相関の絶対値|pij|が最小となる (i, j)」は,次の関数で求まることを確認せよ. ただしグラフGに対応する行列を amat とおく. グラフの辺のうち,偏相関の絶対値が最小となる辺を求める select.ij <- function(P,amat){ p <- nrow(P); minabsP <- Inf for(i in (2:p)){
for(j in (1:(i-1))){
if(amat[i,j] == 1 && abs(P[i,j]) < minabsP){ minabsP <- abs(P[i,j]); i0 <- i; j0 <- j }}} c(i0,j0) } 以下, 例で見てみよう. 標本から計算される相関行列と偏相関行列は表 1 で与えられていたが, 下に再掲する. それぞれ M0, P0 とおく. M0 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.41 0.485 0.711 1 S5 0.389 0.436 0.665 0.607 1 P0 S1 S2 S3 S4 S5 S1 —— S2 0.328 —— S3 0.229 0.282 —— S4 -0.001 0.078 0.432 —— S5 0.026 0.02 0.357 0.253 —— モデル P0 のAICは 0である： AIC0 = 0. 偏相関行列 P0 において絶対値が最小なものは, (4,1)の−0.001である. なので辺(4,1)を除い たグラフを検討する. fitConGraph で推定された相関行列と,そこから計算した偏相関行列およ びAIC は,
M1 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.41 0.485 0.711 1 S5 0.389 0.436 0.665 0.607 1 P1 S1 S2 S3 S4 S5 S1 —— S2 0.328 —— S3 0.229 0.282 —— S4 0 0.078 0.432 —— S5 0.025 0.02 0.357 0.253 —— AIC1 =−2 となる. AIC1 ≤ AIC0 なので先に進む. 次に P1 の中で絶対値最小の偏相関係数は, (5,2)の0.02である(P0 でなくP1 に注目している ことに注意). そこで, (5,2)を除いたグラフを検討し, M2 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.411 0.485 0.711 1 S5 0.389 0.426 0.665 0.607 1 P2 S1 S2 S3 S4 S5 S1 —— S2 0.329 —— S3 0.225 0.289 —— S4 0 0.082 0.428 —— S5 0.032 0 0.362 0.254 —— AIC2=−3.96 が得られる. AIC2≤ AIC1 なので,さらに次に進む. 以下同様にして, P2 の中では(5, 1) 成分が絶対値最小 (0.032)なので, M3 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.406 0.485 0.711 1 S5 0.368 0.419 0.665 0.607 1 P3 S1 S2 S3 S4 S5 S1 —— S2 0.33 —— S3 0.24 0.286 —— S4 0 0.085 0.424 —— S5 0 0 0.372 0.255 —— AIC3=−5.86 P3 の中では (4, 2)成分が絶対値最小 (0.085)なので, M4 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.546 0.61 1 S4 0.388 0.433 0.711 1 S5 0.363 0.405 0.665 0.607 1 P4 S1 S2 S3 S4 S5 S1 —— S2 0.331 —— S3 0.235 0.327 —— S4 0 0 0.451 —— S5 0 0 0.364 0.256 —— AIC4 =−7.1 P4 の中では (3, 1)成分が絶対値最小 (0.235)なので, M5 S1 S2 S3 S4 S5 S1 1 S2 0.553 1 S3 0.337 0.61 1 S4 0.24 0.433 0.711 1 S5 0.224 0.405 0.665 0.607 1 P5 S1 S2 S3 S4 S5 S1 —— S2 0.465 —— S3 0 0.391 —— S4 0 0 0.464 —— S5 0 0 0.374 0.256 —— AIC5 = 0.208
を得る. AIC5 > AIC4 となったので, P4 を共分散選択の結果とする. このときのグラフが図 2 の ようになる. 問 6. 共分散アルゴリズムを実装し,上記の計算結果を確認せよ．ただし, アルゴリズム全体をい きなり一つの関数にまとめず,最初は P0 からP1 を求める部分だけ実装するなどして,動作を確 認しながら作業を進めること. また，web ページに置いてある付録を読み，「無向独立グラフに対 応する相関行列」とは何か，理解すること． （実習２日目はここまでを目標とする.)
