[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
●平均値(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>>
標準偏差の計算
標本ベクトル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>
変量の標準化(標本ベクトルの標準化)
標準化した標本ベクトル
│ 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となり、
誤差ベクトル>相関ベクトルとなり意味なくなる。
標準化したデータ行列
│ 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
情報群を多変量として整理する。(1)その中で、式による予測および関係式の
発見: 重回帰・正準相関分析。I型。(2)変量をもとに基本式をつくる:
判別分析。II型。(3)基準となる変量をさがす: 主成分分析。III型。がある。
いままで、主として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をつかった統計的分析を多変量分析という。
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
となる。
合成変量ベクトル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
●ベクトルの線形変換 平面上のベクトル
│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が
旧い座標で見たらどうなるか。
となる。固有値は、変換した新しい各座標軸の大きさの係数。
●(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と同じ。)
関係式の発見につかう。説明変数
[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 の推定値とする。
これが最小二乗法である。
●第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 │
●分散分析
目的変量(測定値)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
└──────────────────────────────┘
●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次式によるあてはめ
m+n種目の変量
[X1 X2・・・Xm Y1 Y2・・・Yn]
があるとき、それを2群にわけ、2つの合成変量
_
X=X・α=α1X1+α2X2+・・・+αmXm
_
Y=Y・β=β1Y1+β2Y2+・・・+βnYn
をつくる。X,Y 間の相関が最大となるように、係数 αi、βj を決めのが、正準
相関分析である。
[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 |
└─────────────────┘
m項目からなる多変量
A=[A1 A2 ・・・・・Am]
がある。これを
Ai=[A1i A2i ・・・ Ami] (i=1, k)
のようにi群に分ける。このとき合成変量(線形判別関数)
Zi=Ai*β=β1Ai+・・・・・・+βmAm (但しββ'=1)
を作り、群間でZiが最大に異なるようにβを決める。
今、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 │
└───────────────┘
となる。
2群の中央値
_ 1 _ _ 1 _ _
Z= ── (Z1+Z2)= ── β'(A1+A2)
2 2
となる。一つの観測値は、2群の一方から得られる確率が等しいから。判別関数と
中央値の差により個別データの当てはまりの良い/悪いがわかる。すなわち、
┌─────────────────┐
| _ |
| β'X − Z≧0 → グループ I |
| _ |
| β'X − Z<0 → グループ II |
└─────────────────┘
である。
●マハラノビス(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 │
└──────────────────────────────┘
変量の総合・整理に役立つ。変量ベクトル(行列)
│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│
を求めることは容易。(注終り)
----------------------------------------------
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
が最大となる。
●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, |
└─────────────────────────────┘
m列 →
D(0,0) ........ D(0,m)
. .
n↓ . .
行 . .
D(n,0) ........ D(n,m)
_corr_mat.........corr
_jacobi...........jac
_mean_sgm.........mean_sgm
_normalize........norm
●二項分布 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
−∞
となる。