IX. 統計量のベクトル的見方

9−1。ベクトルの射影

[X]の 単位ベクトル[a]への射影は [Xa]=Πa[X] である(図9−1)。ただし、 [a][a]' Πa= ──── [a]'[a] これを、射影子とよぶ。 <図9−1> ‖Xa‖2 = Xa'・Xa = Xa・Πa・X [X]の[a]、および直交している[a]⊥への射影: [Xa⊥] = [ [I]−[Πa] X ただし、 │ 1 │ I= │ 1 │ │ ・ │ │ 1 │ これを、conベクトル(行列)という。 <Prog 9-1> ベクトルXのベクトルaへの射影 ┌───────────────────────┐ │Project vector: to vector: makes vector: │ │ 1.000000 2.000000 0.666667 │ │ 1.000000 2.000000 0.666667 │ │ 0.000000 2.000000 0.666667 │ └───────────────────────┘ X は、 X =Xa +Xa⊥ である。ただし、Xa⊥は、Xaと直交するベクトル。[a]'を掛けると、 [a]' X =[a]'Xa +[a]' Xa⊥ ‖ 0 である。一方、Xa は [a]上にあるから、 Xa= (α)[a] と書ける。[a]'を掛けると、 [a]' Xa =(α)a'a 従って [a]' X =(α)a'a それで a' α= ──── X a'a a a' Xa = [────] X=ΠaX a'a

9−2。平均値と標準偏差

●平均値(1変量) 変量=平均値+バラツキ 例 国語の得点=[難しさ]+[一生徒のバラツキ] 変量のとった値全体として[標本ベクトル] │ a1 │ │ μ │ │ b1 │ │ a2 │ │ μ │ │ b2 │ │ ・ │ = │ ・ │ + │ ・ │ │ ・ │ │ ・ │ │ ・ │ │ an │ │ μ │ │ bn │ 又は A U B (標本ベクトル)=(平均ベクトル)+(偏差ベクトル) ここで平均ベクトルを求めるには、‖B‖(偏差)を最小とする μを求めればよい。これはAのconベクトルI への射影となる。すなわち、 AI=ΠI・A I・I' 1 = ─── ・ A = I ・(───)・(I'・A) I'・I I'・I 1 = I・ ──・(a1+・・・+ai+・・・+an)} n そして、conベクトルの何倍かが平均値となる。 ●標準偏差 standard deviation σ │ b1 │ │ a1−μ │ │ b2 │ │ a2−μ │ 偏差 B=│ ・ │=標本ベクトル A−平均ベクトル U= │ ・ │ │ ・ │ │ ・ │ │ bn │ │ an−μ │ 標準偏差とは、「B の絶対値は、conベクトルの何倍か」、すなわち ‖B‖=σ‖I‖ とおいたときの σ。 --------------------------------------- ‖B‖2 1 1 σ2= ──── = ──‖B‖2 = ── (b12+b22+・・・・・・bn2) ‖I‖2 n n 1 = ──{(a1−μ)2+(a2−μ)2+・・・・(an−μ)2} n また、 1 1 σ2= ──‖B‖2 = ── (‖A‖2−‖μ・I‖2) n n 1 = ──‖A‖2−μ2 n と計算出来る。 <Prog 9-2>> 標準偏差の計算

9−3。標本ベクトルの標準化

標本ベクトルB を、平均0、標準偏差1のベクトルA に変換する。 1 │ a1 │ A= ──(B−μ)= │ a2 │ σ │ ・ │ │ an │ ただし、μ$σはBの平均・標準偏差。このとき I'A= a1 +a2 +a3 +・・・・an =0 A'A= a12+a22+a32+・・・・an2=n となる。 <Prog 9-2> 変量の標準化(標本ベクトルの標準化)

9−4。2変量と相関係数

標準化した標本ベクトル │ a1 │ │ b1 │ 但し A'・A=B'・B=n │ a2 │ │ b2 │ I'・A=B=0 A= │ ・ │ B=│ ・ │ │ ・ │ │ ・ │ │ an │ │ bn │ いま、Aを基準として、BがAの何倍 (k) かが、相関係数である。 すなわち、 A−B・k=0 になるように k をきめる。この誤差を最小とすることは、BのA上への射影BA が、Aの何倍(k)かを調べること。したがって、 A・A' A'・B 1 ΠA・B = ──── B = A・ ──── = A・( ── A'・B) A'・A A'・A n 1 k = ── ・A'・B n -------------------------------------------------------------- 以下、A,B 標準化しない前の標本ベクトルとすると、 1 1 1 = ── ・[ ── (A−μA)][ ── (B−μB)] n σA σB 1 (A−μA )' (B−μB ) = ───── ・ ──────────── σA・σB n Σ(ai・bi/n − μA・μB) = ─────────────── σA・σB ------------------------------------------------- 相関係数は標本ベクトルの一次変換で変化せず。 ------------------------------------------------- ここで誤差ベクトルの大きさ ‖B−A・k‖μ2 =‖B‖2 −‖A‖2・k2 = n(1−k2) 1 k < ── = 0.7となる。 2 ‖B−A・k‖2 > ‖A・k‖2となり、 誤差ベクトル>相関ベクトルとなり意味なくなる。

9−5。相関行列:

標準化したデータ行列 │ A1 │ 1 1 │ A2 │ �堯� ── A'A = ── │ ・ │ ・ [A1,A2・・・・An] n n │ ・ │ │ An │ │ A1・A1, A1・A2 ・・・・ A1・An │ 1 │ ・ ・ ・ │ = ── │ ・ ・ ・ │ n │ ・ ・ ・ │ │ An・A1 ・・・・・・・・ An・An │ │ 1 ・・・・・・・・・・・γ1n │ │ γ21・ │ │ ・ │ = │ ・ │ │ ・ │ │ ・ │ │ γn1・・・・・・・・・・・ 1 │ <Prog 9-2> 標準化した相関行列(→多変量): CORR.C

X. 多変量分析の基礎

10−1。多変量分析の大要

情報群を多変量として整理する。(1)その中で、式による予測および関係式の 発見: 重回帰・正準相関分析。I型。(2)変量をもとに基本式をつくる: 判別分析。II型。(3)基準となる変量をさがす: 主成分分析。III型。がある。

10−2。多変量

いままで、主として1種目の標本ベクトル(n個のデータ)を扱ってきた。 │a1.│ │a2.│ A..=│ ・ │ │ ・ │ │an.│ 多変量分析では、m種目の標本ベクトル: │a11│ │a12│ │a1m│ │ ・ │ │ ・ │ │ ・ │ A.1=│ ・ │ A.2=│ ・ │ ・・・ A.m=│ ・ │ │ ・ │ │ ・ │ │ ・ │ │an1│ │an2│ │anm│ それを纏めた「データ行列」 A = [A.1, A.2・・・・・A.m] を扱う。Aをつかった統計的分析を多変量分析という。

10−3。合成変量

n個のm種別変量ベクトル(データ行列) │a1i│ │ ・ │ [A]=[A1, A2 ・・・・・・Am] Ai= │ ・ │ (i=1, m) │ ・ │ │ani│ がある。各種目の重みベクトル │l1 │ │・ │ L =│・ │ │・ │ │lm │ と線形結合させ、 │z1│ │z2│ [Z]=│・ │=[A]・[L]=l1・[A1] +l2・[A2] +・・・lm・[Am] │・ │ │zn│ をつくるとき、Zを合成変量ベクトルという。その平均、分散は E(Z)=I' Z (Iは、conベクトル である) 1 1 σ2(Z)= ── (Z'Z) =L'(── A'A) L=L'ΣA L n n である。 通常、Z は標準化されている事が多い。標準化するには、 Z−I'Z (I−I')A・L X= ────── = ────────── σ √ L' ΣAL ̄ かくして、標準化されたXについて、平均、分散は E(X)=I' X=0 1 σ2(X)= ── (X'X) =1 n となる。

10−4。合成変量の構造

合成変量ベクトルX と、その構成変量ベクトルAi(i=1, m) の集まった行列A との間の相関係数を合成変量の構造という。両者、標準化してものについては、 │γ1 │ │ ・ │ 1 R =│ ・ │= ── A'・X │γm │ n である。 <Prog 10-1> 合成変量とその構造 (2項目、4データ) *** Data 1 Data 2 ************************ 2.800000 5.000000 μ: 1.000000 4.400000 2.45 7.48 5.000000 15.500000 σ: 1.000000 5.000000 1.65 4.64 normalized Data: 0.212708 -0.533439 μ: -0.881219 -0.662757 -0.00 0.00 1.549730 1.729635 σ: -0.881219 -0.533439 1.00 1.00 *** Weight Vectors ************************* L 1 L 2: 3.000000 2.000000 5.000000 1.700000 *** Z ************************************* Z1(L1) Z2(L2) 33.400000 14.100000 μ: 25.000000 9.480000 44.73 17.61 92.500000 36.350000 σ: 28.000000 10.500000 27.75 10.96 normalized Data: -0.408157 -0.320137 μ: -0.710896 -0.741815 -0.00 -0.00 1.721829 1.710669 σ: -0.602775 -0.648717 1.00 1.00 ******** Structures ************************ Z1 Z2 0.934796 0.952085 data 1 0.997140 0.991825 data 2

10−5。行列の固有値及び固有ベクトルの意味

●ベクトルの線形変換 平面上のベクトル │x│ W =│ │ │y│ と対称行列 │ a11 a12 │ A = │ │ │ a21 a22 │ の積 Z =WA のAは、Wの線形変換をあらわす。しかし、それがどんなベクトルに移る(変換) か直観的にぴんとこない。いま、 WoA=Woλo (但し、Wはベクトル。λは実数。) W1A=W1λ1 を満足するλx、Wx があるとする。そのとき、Wは W=aWo+bW1 と書くことができ、それのAによる変換は Z=WA=(aWo+bW1)A = aWoA+bW1A = λo(aWo)+λ1(bW1) のように線形結合式で表わすことが出来る。その変換は、Wo方向にλo倍、W1方向 にλ1倍の和になることがわかる。このとき、λを行列Aの固有値、Wを 固有ベクトルという。 ●行列Aと固有値の関係 二次形式(3次元を使って)をつかって、対応する 事が出来る。一般に │x│ W = │ │ │y│ で、対称行列 │ a11 a12 │ A = │ │ │ a21 a22 │ Z =W'AW =a11x2+2a12xy+a22y2 となり、二次曲面を示す。Z=一定のとき x,y は x'-y'軸 を径とした楕円と なる。 <図10ー1> それで、x'−y'軸にたいしては │x1│ │λ1 0 │ W1=│ │ D =│ │ │y1│ │ 0 λ2│ とすると Z=W1'DW1 =λ1x2+λ2y2 となる。各点のx−y軸、x'−y'軸での関係は W=[P1 P2] W= P・W1 ( 但し、P1'P2=0 ) それで Z=W'PW=P(PW1)'A(PW1)= W1'P'APW1 従って P'AP=D or AP=DP となる。それで、Dの対角要素がAの固有値、Pが正規固有ベクトルとなる。 <図10−2> ●固有値ベクトルの意味 結局、新しい座標で見たW1が、旧い座標で見たら どうなるかという事である。新しい座標で、 │ x1 │ │ 1 │ │ 0 │ W1=│ │ │ │ │ │ │ y2 │ │ 0 │ │ 1 │ の各点は、もとの座標では、 │ x │ │ cosθ │ │ -sinθ │ W= │ │ │ │ │ │ │ y │ │ sinθ │ │ cosθ │ 古い軸から新しい軸への角度をθとすると │ cosθ -sinθ │ P=│ │ │ sinθ cosθ │ を変換行列としたとき、 W = P W1 ↑ ↑ | 新しい座標で見たW1が 旧い座標で見たらどうなるか。 となる。固有値は、変換した新しい各座標軸の大きさの係数。

10−6。行列の固有値及び固有ベクトルの算出

●(a)式を解く │ a11 a1m │ │ ・ │ A = │ ・ │ │ ・ │ │ a1m amm │ のとき AL = λL となるような │ l1 │ │ ・ │ λ, L =│ ・ │ │ ・ │ │ lm │ を求める。まずλをきめる。 AL−λL=0 (A−λI)L=0 L≠0 の条件下で det(A−λI)L=0 これからλをきめる。λがきまったら AL−λL=0 に代入してLをきめる。 <例> │2 1│ │ │ の固有値、固有ベクトルを求める。 │4 -1│ 固有値は、 │2 1│ │1 0│ │2-λ 1 │ det(│ │−λ│ │)=det│ │=0 │4 -1│ │0 1│ │4 -1-λ│ から (2-λ)(-1-λ)-4=0 それで、 λ=3, λ=-2 となる。λ=3 に対する固有ベクトルは、 │2 1│ │l1│ │l1│ │ │*│ │=3 *│ │ │4 -1│ │l1│ │l2│ から l1−l2=0 それで、K=常数とすれば、 │1*K│ L=│ │ │1*K│ となる。通常、一番簡単な K=1 としたものを固有ベクトルとして採用する。 ●(b)Jacobi法 対称行列の対角変換を行うことによって固有値、固有ベクトルを求める方法である。 │1 ・ ・r1m │ │ 1 ・ ・ 0 │ │r21・ ・ │ │ 1 ・ ・ │ │・・・・・・・・・ │ k行 │・・・C・・・-S・・・│ │ ・ ・ │ │ ・ ・ │ Σ=│ ・ ・ │ P= │ 0 ・ ・ 0 │ │ ・ ・ │ │ ・ ・ │ │・・・・・・・・・ │ l行 │ ・・ S・・・ C・・・│ │ ・ ・ │ │ ・ ・ 1 │ │rm1・ ・ 1 │ │ 0 ・ ・ 1 │ k列 l列 k列 l列 rij=rji C=COSθ S=SINθ 固有値:P'ΣPを行ってΣを対角化する。 k,l(又はl,k)部が 0 となる 変換は 2akl tan 2θ= ────── (akk≠all) akk−all cos 2θ= 0 (akk=all) これらの式をみたす。θは、 x=2akl、 y=akk−all と置いたとき、 y≧0 ならば R=√x2+y2 ̄ y<0 ならば R=−√x2+y2 ̄ 1 y cosθ= √──(1 + ──) ̄ 2 R sinθ= √1− cos2θ ̄ となる。この変換を、全要素に行う。そのとき、主対角要素が固有値となる。 固有ベクトル:上記変換をするごとに、単位行列E*Pを行う。これが 固有ベクトルとなる。 <Prog 10-2> 固有値・固有ベクトルをJacobi法で算出するプログラム。 下記のデータを入力すると ┌───────────┐ │3 3 │ │ 1.0 2.0 0.0 │ │ 2.0 2.0 2.0 │ │ 0.0 2.0 3.0 │ └───────────┘ 次の答えが結果する。 ┌─────────────────┐ │Data Matrix: │ │ 1.000000 2.000000 0.000000 │ │ 2.000000 2.000000 2.000000 │ │ 0.000000 2.000000 3.000000 │ │Eigen Values: │ │ 5.000000 2.000000 -1.000000 │ │ │ │Eigen Vectors: │ │ -0.333288 0.666667 -0.666689 │ │ -0.666712 0.333333 0.666621 │ │ -0.666644 -0.666667 -0.333379 │ └─────────────────┘ ●(c)パワー法 最大の固有値・固有ベクトルの算出 すべての正方行列に対応する。m次 行列Xの固有値を λ1, λ2・・・・λm (λ1>λi) とし、それに対応する固有ベクトルを L1,・・・・Lm とする。任意のベクトルは B=X1・L1+X2・L2+・・・・・+Xm・Lm と書ける。両辺にXをかけると X・B=XX1L1+XX2L2+・・・+XXmLm X・Li=λi・Li であるから、 =λ1X1L1+λ2X2L2+・・・+λmXmLm さらにXをかけると X2B=λ12X1L1+λ22X2L2+・・・+λm2・Xm・Lm n回かけると、 D(n)=XnB=λ1nX1L1+λ2nX2L2+・・・+λnm・Xm・Lm λ2 n X2 =λ1nX1{L1+(───) (───)L2+・・・ } λ1 X1 ┌───────────┐ │D(n)=λ1n・X1・L1 │ (係数λ1n・X1とL1'L1=1となるように選ぶ) └───────────┘ D(n)'・D(n)=λ12nX12L1'・L1 D(n)'・D(n−1)=λ12n−1X12L1'・L1 n>>0ならば ┌──────────────┐ │ D(n)'・D(n) │ │λ1= ────────── │ │ D(n)'・D(n−1) │ └──────────────┘ 第2の固有値、固有ベクトル X1=X−λ1・L1・L1' (但しL1'L=1) X1の固有値、固有ベクトルをλ1・Lとする。 第2の固有値 λ2 第2の固有ベクトル L2=(λ1−λ2)L−λ1*(L1'・L)*L1 となる。 <Prog 10-3> パワー法による行列の固有値、固有ベクトルの算出。 (操作は、Prog 10-1と同じ。)

XI. 重回帰分析(multiple regression analysis)

11−1。重回帰分析の大要

関係式の発見につかう。説明変数 [X1, X2,・・・・Xp] で目的変数Yを説明しようとして次の模型 Y=βo+β1X1+β2X2+・・・・・βpXp+ε を考える。Xiは区間 [ci, di] 内での1つの母集団とする。εは誤差といわれ、 (1)E(ε)=0、 Var(ε)=σ2の分布 (2)εの分布は正規分布 を仮定。このとき E(Y)=βo+β1X1+β2X2+・・・・・・・βpXp Var(ε)=σ2 (Xi=xiをとる時) ----------------------- 実際に得られたn組のデータ(大きさnの標本)の内、k番目のデータは yk= bo+b1xk1+b2xk2+・・・・・bpxkp+ek (k=1,・・n) という関係が成立する。このとき、誤差項ekの平方和 n Σek2=Σ{yk−(bo+b1xk1+・・・・・・・bpxkp)}2 k=1 を最少になるような係数bo, b1・・・bp を βo,β1・・・βp の推定値とする。 これが最小二乗法である。

11−2。回帰式への当てはめ

●第1の解法 │ y1 │ │1 x11 x21・・・xp1 │ │ bo │ │ y2 │ │1 x12 x12 │ │ b1 │ Y=│ ・ │ X=│ ・ │ b=│ ・ │ │ ・ │ │ ・ │ │ ・ │ │ yn │ │1 x1n x2n・・・xpm │ │ bp │ とすれば S=Σek2 =(Y−Xb)'(Y−Xb) =YY'−2b'X'Y+b'X'Xb 最小の条件は ∂S ── =−2X'Y+2X'Xb=Zer ∂b 従って X'Xb=X'Y これを正規化方程式(normal equation)という。 ┌───────────┐ │ b=(X'X)−1 X'Y │ └───────────┘ ●第2の解法 また偏差積和行列をつかっても求められる。 │a11・・・a1p │ │ │ _ _ A =│ Σ │ =(X−X)'(X−X) │ │ │ap1・・・app │ │a1y│ │ │ _ _ Ay =│ │=(X−X)'(Y−Y) │ │ │apy│ _ _ _ │b1│ │x11 ・x1p│ │x1 ・xp│ │y│ │y│ │ ・│ │ ・ │ _ │ ・ │ │・│ _ │・│ b=│ ・│ X=│ ・ │ X=│ ・ │ Y=│・│ Y=│・│ │ ・│ │ ・ │ │_ _ │ │・│ │_│ │bp│ │xn1 ・xnp│ │x1 ・xp│ │y│ │y│ X、Y間に一次式の関係があれば、 Ay=A・b │A│≠0 となる。なんとなれば、 0 _ _ ‖ (Y−Y)=(X−X)b +[ek] _ _ _ _ (X−X)'(Y−Y) = (X−X)'(X−X)b Ay A それで b=A−1Ay _ _ _ _ _ bo=y−(b1x1+b2x2・・・・・・+bpxp) ●第3の解法 下の式をとく。 Y=b2X2+b1X1+bo │ΣX22 ΣX2X1 ΣX2 │ │ b2 │ │ ΣX2Y │ │ │ │ │ │ │ │ΣX1X2 ΣX12 ΣX1 │ *│ b1 │ = │ ΣX1Y │ │ │ │ │ │ │ │ΣX2 ΣX1 Σ1 │ │ bo │ │ ΣY │

11−3。重回帰式の有意性の検定

●分散分析 目的変量(測定値)yoの全変動ayy= 回帰推定値(計算値)yの変動+残差eの変動 n _ ayy =Σ(yoi−y)2 i=1 _ =Σ(yoi−yi+yi−y)2 _ _ _ =Σ(yoi−yi)2+Σ(yi−y)2+2Σ(yoi−yi)(yi−y) ↓ ↓ ‖ 残差 推定値の変動 Σe=0 あるいは、 _ _ ayy = (Yo−Y)'(Yo−Y) _ _ = (Yo−Y)'(Yo−Y)+(Y−Y)'(Y−Y) ↓ ↓ 残差の変動(EV) 回帰推定値の変動(RV=Ay'・b) RV >> EV ならば有効。 ------------------------------------------------- │yo1 │ │y1 │ _ Yo= │ ・ │ Y=│・ │ Y=Yo'Yo │ ・ │ │・ │ │yon │ │yn │ 測定値 計算値 <図11−1> 推定値・残差・平均・全変動・回帰推定値の変動の説明図 ┌──────┬───────────┬───┬─────────┐ │ 変動要因 │ 平方和 │自由度│ 不偏分数 │ ├──────┼───────────┼───┼─────────┤ │ │ n _ │ │ │ │ 全変動 │ayy=Σ(yoi−y)2 │ N−1 │ │ │ │ i=1 │ │ │ ├──────┼───────────┼───┼─────────┤ │ │ n _ │ │ │ │ │ RV=Σ(yi−y)2 │ │ │ │ 回帰推定値 │ i=1 │ p │ VR= RV/p │ │ の変動 │ p │ │ │ │ │ =Σakybk │ │ │ │ │ k=1 │ │ │ ├──────┼───────────┼───┼─────────┤ │ │ n │ │ EV │ │ │ EV=Σ(yoi−yi)2 │ │ Ve= ──── │ │ 残差の変動 │ i=1 │ │ n−p−1 │ │ │ =ayy−RV │ │ │ └──────┴───────────┴───┴─────────┘ 分散比 F=VR/Ve は、自由度(p, n−p−1) のF−分布をする。 Ho: 求めた重回帰式はyの推定に何ら役立たない。 H1: accept F → F(p, n−p−1;α)ならば、危険率 α%でHoを棄却。 ●重相関係数 目的変量(標本値、測定値)yoi ←→ 回帰推定値 yi n _ _ Σ(yoi−y)(yi−y) i=1 R= ──────────────── n _ n _ √( Σ(yoi−y)2 Σ(yi−y)2 ) ̄ i=1 i=1 寄当率(proportion) RV(回帰による変動) = ──────────── ayy(全変動) <Prog 11-1> 重相関式へのあてはめと寄与率の算出:do_regr ┌───────────┐ │9 3 │←9データ、3=(1目的変数+2説明変数) │1.25 16. 96.5 │ │3.3 2.4 43.0 │ │1.2 4.6 33.5 │ │2.1 6.0 49.9 │ │3.4 7.2 70.3 │ │1.5 11.0 70.7 │ │4.0 14.0 103.0 │ │3.0 12.7 95.3 │ │2.5 17.0 120.6 │ │ │ │2 │← 0,1,2 の内、2列目が目的変数 └───────────┘ ↓ ┌──────────────────────────────┐ │X: Y: Y-expected: │ │ 1.000000 1.250000 16.000000 96.500000 97.819062 │ │ 1.000000 3.300000 2.400000 43.000000 40.790349 │ │ 1.000000 1.200000 4.600000 33.500000 35.204854 │ │ 1.000000 2.100000 6.000000 49.900000 50.380658 │ │ 1.000000 3.400000 7.200000 70.300000 67.815433 │ │ 1.000000 1.500000 11.000000 70.700000 72.634183 │ │ 1.000000 4.000000 14.000000 103.000000 109.939648 │ │ 1.000000 3.000000 12.700000 95.300000 94.471888 │ │ 1.000000 2.500000 17.000000 120.600000 113.743924 │ │ │ │Regr Parameters(b0,b1,.....): │ │ 0.058091 8.375297 5.455741 ←(Y=bo・X1+b1・X2+b2) │ │ │Corr: 0.991886, Proportion: 0.983838 │ │F: 101.454082 (3, 5) ←p(3), n(9)-p(3)-1 └──────────────────────────────┘

11−4。重回帰式への当てはめ(続)

●2変数1次式によるあてはめ(#3の解法) y=b2X2+b1X1+bo │X22 X2X1 X2 │ │ b2 │ │ΣX2Y │ │ │ │ │ │ │ │X1X2 X12 X1 │ * │ b1 │ = │ΣX1Y │ │ │ │ │ │ │ │X2 X1 1 │ │ bo │ │ΣY │ <Prog 11-2> 2変数1次式によるあてはめ ●1変数3次式によるあてはめ Z=a3X3+a2X2+a1X+ao │ΣX6 ΣX5 ΣX4 ΣX3│ │a3 │ │ΣX3・Z│ │ │ │ │ │ │ │ΣX5 ΣX4 ΣX3 ΣX2│ │a2 │ │ΣX2・Z│ │ │*│ │ = │ │ │ΣX4 ΣX3 ΣX2 ΣX │ │a1 │ │ΣX・Z │ │ │ │ │ │ │ │ΣX3 ΣX2 ΣX Σ1 │ │ao │ │ ΣZ │ <Prog 11-3> 1変数3次式によるあてはめ ●2変数 3次式によるあてはめ Z=a3X3+a2X2+a1X+b3y3+b2y2+b1y+b0 を最小二乗法で a3, a2, a1, b3, b2, b1, b0を求める。 │ΣX6 ΣX5 ΣX4 ΣX3Y3 ΣX3Y2 ΣX3Y ΣX3│ │ΣX5 ΣX4 ΣX3 ΣX2Y3 ΣX2Y2 ΣX2Y ΣX2│ │ΣX4 ΣX3 ΣX2 ΣXY3 ΣXY2 ΣXY ΣX │ │ΣX3Y3 ΣX2Y3 ΣXY3 ΣY6 ΣY5 ΣY4 ΣY3│ │ΣX3Y2 ΣX2Y2 ΣXY2 ΣY5 ΣY4 ΣY3 ΣY2│ │ΣX3Y ΣX2Y ΣXY ΣY4 ΣY3 ΣY2 ΣX │ │ΣX3 ΣX2 ΣX ΣY3 ΣY2 ΣY Σ1 │ │ a3 │ │ ΣX3・Z │ │ a2 │ │ ΣX2・Z │ │ a1 │ │ ΣX ・Z │ *│ b3 │= │ ΣY3・Z │ │ b2 │ │ ΣY2・Z │ │ b1 │ │ ΣX ・Y │ │ bo │ │ ΣZ │ E*E1=E0 <Prog 11-4> 重相関分析の例: 2変数 3次式によるあてはめ

XII. 正準相関分析 Canonical correlation

12−1。分析の大要

m+n種目の変量 [X1 X2・・・Xm Y1 Y2・・・Yn] があるとき、それを2群にわけ、2つの合成変量 _ X=X・α=α1X1+α2X2+・・・+αmXm _ Y=Y・β=β1Y1+β2Y2+・・・+βnYn をつくる。X,Y 間の相関が最大となるように、係数 αi、βj を決めのが、正準 相関分析である。

12−2。分析の操作

[X, Y] の相関行列を X Y ↓ ↓ │ Σ11 Σ12 │ ←X Σ=│ │ │ Σ21 Σ22 │ ←Y であらわし、相関係数を COV(X,Y) ρc= ─────── ��V(X)・V(Y) であらわす。但し、 V(X)=α'Σ11α=1 (標準化しているから) V(Y)=β'Σ22β=1 _ _ COV(X,Y)=COV(Xα,Yβ)=α'Σ12β である。そのとき、 ρc=α'Σ12β となり、正準相関分析は、これを最大化することと等価となる。最大になるのは、 Σ12β−λΣ11α=0 Σ12α−λΣ12β=0 のときである。(證明略)。これを解くと、 (Σ11-1・Σ12・Σ22-1・Σ21−λ2I)α=0 (Σ22-1・Σ21・Σ11-1・Σ12−λ2I)β=0 第二式についてλ2・βを求める。これは�� Σ22-1・Σ21・Σ11-1・Σ12 の固有値λ2、固有ベクトルβ を求めること。である。つぎに、 �熙�11-1Σ12β α= ──────── λ でαを求める。またそのときの相関形数 ρc=α'Σ12β=λα'Σ11α=λ ( α'Σ11α=1 であるから ) である。 <Prog 12-1> 正準相関分析 <相関分析の例> 相関行列を入力: ┌───────────────────────────┐ |5 5 | | 1.000000 -0.249030 -0.210051 -0.401331 0.361817| | -0.249030 1.000000 0.288697 0.045946 -0.217097| | -0.210051 0.288697 1.000000 0.652113 -0.902671| | -0.401331 0.045946 0.652113 1.000000 -0.881982| | 0.361817 -0.217097 -0.902671 -0.881982 1.000000| |2 |←X,Y境界 └───────────────────────────┘ ↓ ┌─────────────────┐ |select boundary of X,Y : | | 1.000000 -0.249030 | | -0.249030 1.000000 | | | | -0.210051 -0.401331 0.361817 | | 0.288697 0.045946 -0.217097 | | | | 1.000000 0.652113 -0.902671 | | 0.652113 1.000000 -0.881982 | | -0.902671 -0.881982 1.000000 | | | | -0.210051 0.288697 | | -0.401331 0.045946 | | 0.361817 -0.217097 | | | |λmax of Y: 0.466126 | |Y: | | -0.215906 -0.835518 0.505267 |βoyo+β1y1+β1y2=αoxo+α1x1 |X: |が関係式式 | -0.060655 1.259028 | └─────────────────┘

XIII. 判別分析(discriminant analysis)

13−1。分析の大要

m項目からなる多変量 A=[A1 A2 ・・・・・Am] がある。これを Ai=[A1i A2i ・・・ Ami] (i=1, k) のようにi群に分ける。このとき合成変量(線形判別関数) Zi=Ai*β=β1Ai+・・・・・・+βmAm (但しββ'=1) を作り、群間でZiが最大に異なるようにβを決める。

13−2。分析の操作

今、2群(i=2)とする。 A1=[A11 A12 ・・・ A1m] A2=[A21 A22 ・・・ A2m] 平均 _ _ _ _ A1=[A11 A12 ・・・ A1m] _ _ _ _ A2=[A21 A22 ・・・ A2m] このとき、群間変動/群内変動 を最大とするようなβを求める。 ●群間変動(各グループ内の得点の平均の間の変動) _ _ (A1− A2)2 _ _ =(β'A1−β'A2)2 _ _ _ _ =β'(A1−A2)(A1−A2)'β ●群内変動 2 _ _ Σ(Ai−Ai)(Ai−Ai)' i=1 _ _ =β'Σ1β+β'Σ2β 但し Σx=(Ax−Ax)'(Ax−Ax) =β'(Σ1+Σ2)β 1 S= ───── (Σ1+Σ2) n1+n2−2 すると、 群内変動=β'[(n1+n2−2)S]β ●群間変動/群内変動 _ _ _ _ 群間変動 1 β'(A1−A2)(A1−A2)'β λ= ────── = ──────・─────────────── 群内変動 n1+n2−2 β'Sβ ●最大化は、∂λ/∂β=0 すなわち _ _ C(A1−A2)−Sβ=0 C=const ●結局 ┌───────────────┐ │ _ _ │ │β=C・S-1(A1−A2) 但し C=1 │ └───────────────┘ となる。

13−3。結果による個々データの分類

 2群の中央値 _ 1 _ _ 1 _ _ Z= ── (Z1+Z2)= ── β'(A1+A2) 2 2 となる。一つの観測値は、2群の一方から得られる確率が等しいから。判別関数と 中央値の差により個別データの当てはまりの良い/悪いがわかる。すなわち、 ┌─────────────────┐ | _ | | β'X − Z≧0 → グループ I | | _ | | β'X − Z<0 → グループ II | └─────────────────┘ である。

13−4。判別関数の有用性の検定

●マハラノビス(Mahalanobis)のD2(またはGeneraliged distance) ┌──────────────────┐ │ _ _ _ _ │ │ D2=(A1−A2)'・S-1(A1−A2) │ └──────────────────┘ ●判別関数に関する汎距離 Do2=(C'ΔX)'(C'SC)-1(C'ΔX) _ _ ΔX=X1−X2 C'=[c1 c2・・・cp] ●有効性の判定。 m=n0・n1/((n0+n1)(n1+n2−2)) とするとき、 n0+n1−p−1 m(D2-Do2) F= ───────── ・ ──────── p−1 1+mDo2 >Fαi;p-1,n1+n2-p-1 ならば、あわないとして棄却する。 <Prog13-1> 判別分析 <判別分析の例> 項目1、項目2、分類番号からなるデータ ┌────────────────────┐ │23 3 │ │750 590 1 670 600 0 │ │360 600 1 560 550 0 │ │720 750 1 540 520 0 │ │540 710 1 590 540 0 │ │570 700 1 590 700 0 │ │520 670 1 470 600 0 │ │590 790 1 560 540 0 │ │670 700 1 540 630 0 │ │620 730 1 500 600 0 │ │690 840 1 │ │610 680 1 │ │550 730 1 │ │590 750 1 │ └────────────────────┘ ↓ ┌──────────────────────────────┐ │μ0: 598.46 710.77 │ │ │ │μ1: 576.00 596.00 │ │ 0 -2.535992 0 │ │ 1 -0.562820 0 │ │ 2 2.220990 1 │ │ 3 1.841339 1 │ │ 4 1.422588 1 │ │ 5 0.770834 1 │ │ 6 3.939170 1 │ │ 7 0.990805 1 │ │ 8 2.074343 1 │ │ 9 4.953463 1 │ │ 10 0.671444 1 │ │ 11 2.376591 1 │ │ 12 2.782309 1 │ │ │ │ 0 0.110125 1 │ │ 1 -1.901349 0 │ │ 2 -2.872464 0 │ │ 3 -3.653754 0 │ │ 4 -3.291215 0 │ │ 5 1.336232 1 │ │ 6 -1.037782 0 │ │ 7 -3.161680 0 │ │ 8 -0.472384 0 │ │ 9 -1.167317 0 │ │Summary: │ │ CORRECT INCORRECT │ │GROUP I : 11 2 │ │GROUP II: 8 2 │ │ │ │β 0: -0.0043178 │ │β 1: 0.0289215 │ │ │ │Zbar: 16.3613216 │ │T=SQR.ACTUAL= 18.2130990 T=SQR.CRITICAL= 12.2850000 │ └──────────────────────────────┘

XIV. 主成分分析 Princeple compoment analysis

14−1。主成分分析の意味

変量の総合・整理に役立つ。変量ベクトル(行列) │a11│ │am1│ A=[A1, A2・・・・・Am] A1=│ ・ │ ・ ・ Am=│ ・ │ │a1n│ │amn│ がある。(但し標準化され、A'A=n)。 │ l1 │ │ ・ │ L= │ ・ │ (L'L=1) │ ・ │ │ lm │ 線形結合し、合成変量ベクトル Z=L*A=l1・A1+l2・A2+ ・・・・・・・・・ lm・Am を作る。(但しL'・L=1)。この時Zの分散 V=Var(Z)=Var(L'*A)=L'ΣL 主成分分析では、V が最大となるLを求める。Σは、Aの分散共分数行列(相関行列) である。 V=L'ΣL=L'λL=λL'・L L'L=1 であるから V=λ(=Σの固有値) となる。それで、このことは、Σの固有値を求めるのと等価になる。Lはそれに 対応する固有ベクトルとなる。 ------------------------------------------- (注1)ラグランジュの未定乗数法 L'L−1=0 (L'L=1) の条件下で V=L'ΣL の最大値を求めるには、 h=L'ΣL+λ(L'L−1) を作り、 dh ── = 2ΣL−2λL=0 dL と同値である。従って L'・L=1 の条件下で ΣL=λ・L となるような λ、Lを求めることになる。 (注2)二次形式 F(X)=X'AX (Aは対物行列) について、1)‖X‖=1のとき、F(X)の最大値、最小値はそれぞれAの固有値の 最大値、最小値と一致する。2)X≠0のとき 1 ──── X'AX ‖X‖2 の最大値、最小値はそれぞれAの固有値の最大値、最小値と一致する。 1)Aの固有値を λ1≦λ2≦・・・・・≦λn とする。適当な直交行列Lによる主軸変換により _ _ _ _ F(X)=λ1 X12 +・・・・+λn Xn'2 (X=L・X) となる。 _ ‖X‖2 =‖X‖2 =1 であるから、 _ _ _ _ _ λ1≦λ1‖X‖2=λ1ΣXj2≦ΣλjXj2 ≦λnΣXj2 =λ1‖X‖2 =λn ∴λ1≦F(X)≦λn 等号はそれぞれ _ _ X'=e1; X'= en のとき成立する。 2) 1 1 1 ───X'AX =(─── X)'A(─── X) ‖X‖2 ‖X‖ ‖X‖ で, 1 ‖─── X‖=1 ‖X‖ であるから(1)と同様である。この方法は、微積分の知識は全然不用。また証明 からわかるようにX'AXが最大(小)となるような │X1│ X=│・ │ │Xm│ を求めることは容易。(注終り) ----------------------------------------------

14−2。主成分分析の別の見方

 n個、m項目のデータ │aio│ │ai1│ A=[Ao A1 ・・・・・Am] Ai=│・ │ (i=0, m) │・ │ │ain│ があり、重みベクトル │lo │ │l1 │ L=│・ │ │・ │ │lm │ をつかって、n個の合成変量データ │zo│ │a1o・lo+・・・amo・lm│ │z1│ │a11・lo+・・・am1・lm│ Zk=│・ │=AL=│ │ │・ │ │ │ │zn│ │a1n・lo+・・・amn・lm│ をつくる。Z とAi との相関係数ri とすれば、 Ao=Z・ro, A1=Z・r2, ・・・・ Am=Z・rm である。このとき、(上述の)「Zの分散を最大にすること」は、 Σ‖Z・ri‖2 ────── Σ‖Ai‖ を最大とする;すなわち、ZとA1, A2・・・ Am とを、出来るだけ高い相関をもつ ようにL をきめるのと等価になる。  それは、また、ZとAiとの誤差ベクトル Σ‖Ai−Z・ri‖2 を最小とすることと等価である。 <図14-1a> 前と対応させると、 1 Zの分散= ──‖Z‖2 − μ2 n が最大となる。

14−3。主成分分析の評価

●Size Factor と Shape Factor  m種目のデータから、固有値、固有ベクトルはm対ある。 Z1= A・L1 Z2= A・L2 ・ ・ ・ ・ Zm= A・Lm が合成変量となる。Z1は、測定値の総合的大きさの順となる(size factor)。これに 反し、Z2は、+−が入ってくる(Shape factor)。個々のデータについての、合成 変量値(得点;score) Z1をX軸にZ2をY軸にプロットすれば、広がりを視覚的に に見ることが出来る。 ●因子負荷量(構造ベクトル) fuctor loading Zと各成分Ai(i=1,・・m)との相関。 │ r1 │ │ r2 │ ΣL R=│ ・ │ = ────── │ ・ │ √L'ΣL ̄ │ rm │ である。ここで、 Σ=λα √L'ΣL ̄=√Σλα ̄ であるから、 R=√λα・L ̄ となる。 ●合成変量の寄与率 Proportion Zkが如何に Ai(i=1,...m) の代表になっているか。AiのZkとの相関。 Zkに対応する固有値λkとする。 <図14-1b> Σ‖Ai‖= n・m i Σ‖Zk・ri‖= n・λk i Σ‖Ai‖の全長中      Σ‖Zk・ri‖   n・λk λk 寄与率= ──────── = ───── ← ───────────      Σ‖Ai‖   n・m λ1+λ2+・・・・λm -------------------------- (注) λ1+λ2+・・・・+λm = tr(A) である。相関行列では、tr(A)=m である。∴Σλi=mとなる。 │a11−X a12 ・・・・・ a1m │ │ ・ ・ │ │ a22−X ・ │ = 0 │ ・ ・ │ │am1 ・ amm−X│ Xn−(a11+a22+・・+amm)Xm−1・・+|A|=0 一方この固有値をλ1,λ2,λ3,・・・λmとすると (X−λ1)(X−λ2)・・・(X−λm)=0 Xn−(λ1+・・・+λm)Xm−1+・・・=0 比較すると λ1+・・・・λm=a11+a22・・・・amm (注終り) ------------------------- <Prog 14-1> 主成分分析 <主成分分析の例> ┌────────────────────────────┐ |5 5 | | 1.000000 -0.249030 -0.210051 -0.401331 0.361817 | | -0.249030 1.000000 0.288697 0.045946 -0.217097 | | -0.210051 0.288697 1.000000 0.652113 -0.902671 | | -0.401331 0.045946 0.652113 1.000000 -0.881982 | | 0.361817 -0.217097 -0.902671 -0.881982 1.000000 | └────────────────────────────┘  ↓ ┌─────────────────────────────┐ |λ: | | 2.882270 1.031113 0.814353 0.251058 0.021207 | |eigen_vect: | | -0.304890 -0.426522 -0.802222 -0.284349 -0.026580 | | 0.198493 0.837153 -0.417826 -0.291373 0.017253 | | 0.519254 -0.105626 -0.373202 0.612240 0.452880 | | 0.518634 -0.278433 0.177795 -0.677747 0.403165 | | -0.573618 0.169033 0.104736 -0.008255 0.794577 | | | | | |lambda; Proportion | | Eigen_vect: | | Factor_loading: | | | |2.882270, 0.576454 | | -0.304890, 0.198493, 0.519254, 0.518634, -0.573618, | | -0.517620, 0.336986, 0.881550, 0.880498, -0.973845, | | | |1.031113, 0.206223 | | -0.426522, 0.837153, -0.105626, -0.278433, 0.169033, | | -0.433106, 0.850077, -0.107256, -0.282731, 0.171642, | | | |0.814353, 0.162871 | | -0.802222, -0.417826, -0.373202, 0.177795, 0.104736, | | -0.723938, -0.377053, -0.336784, 0.160445, 0.094515, | | | |0.251058, 0.050212 | | -0.284349, -0.291373, 0.612240, -0.677747, -0.008255, | | -0.142475, -0.145995, 0.306767, -0.339590, -0.004136, | | | |0.021207, 0.004241 | | -0.026580, 0.017253, 0.452880, 0.403165, 0.794577, | | -0.003871, 0.002512, 0.065951, 0.058711, 0.115711, | └─────────────────────────────┘

XIV. 資料と付録

15−1。データ行列・フォーマット

m列 → D(0,0) ........ D(0,m) . . n↓ . . 行 . . D(n,0) ........ D(n,m)

15−2。統計ライブラリ STAT.LIB

_corr_mat.........corr _jacobi...........jac _mean_sgm.........mean_sgm _normalize........norm

15−3。マトリックス演算ライブラリMAT.LIB

15−4。二項分布と正規分布

●二項分布 B(n,p) において、p=一定で、n → ∞ とすると、分布は連続 となり、正規分布をしめす。 (証明) 成功確率pであるベルヌーイの試行をn回行って、r回の成功の起る確率 B(n,r)=nCr pr(1−p)n-r Rは離散的な偶然変数:r=0, 1, 2 ...n の値をとる二項分布を持つ。 このとき, E[R]=np E[(R−np)2]=npq (但し q=1−p) ◎Rを正規化する (X の平均値は0。分散は1 にする)。 R−np X=──── or R=np−X√npq ̄ √npq ̄ ◎x およびaをある2個の固定された数とすれば Pr[a≦ X ≦x]=ΣB(n,r) である。但し和は不等式 np+a√npq ̄≦ r ≦np+x √npq ̄ の範囲のrのすべての値に対して行う。 ◎n!=nne-n√2πn ̄(1+θn) limθn=0 (スターリング公式) n→∽ を用い、いくらか計算すると、 1 1 log B(n,r)=(np+d+ ──)log n +(nq−d+ ──) logz(n) 2 2 1 1 −(np+d+ ──)log(np+d)−(nq−d+ ──) log(nq−d) 2 2 1 1 +(np+d+ ──)log p+(nq−d+ ──) log(q) 2 2 −log√2πnpq ̄+θn' とかける。ただし d=r−np である。この式は 1 d 1 d log B(n,r)=(np+d+──) log(1+──)+(nq−d+──) log(1− ──) 2 np 2 nq −log√2πnpq ̄+θn' とかける。 ◎ξをαとx との間のある数、すなわちα≦ξ≦x とすれば r=np+ξ√npq ̄ だから d r−np ξ√npq ̄ ── = ──── = ───── np np np である。それで | d | M q |──|≦ ──── √── ̄ , M=max (|α|,|x |) | np | √n ̄ p となる。ここで p と x を固定して、n → ∞ とすればこの不等式の右辺は 0に収束する。それで、 d |──| →0 np それで、上式に z2 log (1+z)= z − ── + ...... |z|<1 (マクローリン展開) 2 を適用できる。 ◎すると、 1 d2 log B(n,r)= − ── ─── − log √2πnpq ̄+θn'' 2 npq となる。すなわち、 1 1 (r−np)2 Prob.[α≦X≦x]=Σ────── exp[− ── ─────](1+θn''') √2πnpq ̄ 2 npq である。 ◎ R−np X= ───── √npq ̄ から、 r−np (r+1)−np xr= ───── , xr+1= ─────── √npq ̄ √npq ̄ および 1 Δxr =xr+1−xr= ───── √npq ̄ とすれば、上記の和は 1 xr^2 Prob.[α≦X≦x]=──── Σ exp(− ────) Δxr (1+θn''') √2π ̄ 2 となる。 ◎この和は積分 1 x x^2 ───── ∫ exp(− ────) dx √2π ̄ a 2 に対する近似和とみなすことができる。 ◎それで、ベルヌーイ分布の極限の分布関数は、 1 x x^2 P(x) = Prob.[X≦x]= ───── ∫exp(────) dX √2π ̄ −∞ 2 となる。これが正規分布(normal distribution)あるいはガウス分布 (Gausssian distribution)のそれに外ならない。 (証明終り)。  対応する度数関数p(X)は 1 x^2 p(X)= ───── exp(────) dX √2π ̄ 2 となり、また、 ∞ ∫ p(X)dX=1 −∞ ∞ E[X]=∫ Xp(X)dX=0 −∞ ∞ E[X2]=∫ X2p(X)dX=1 −∞ となる。