« 最近考えていること(場の理論等覚え書き) | トップページ | 確率と分布関数(11)(区間推定)(終了) »

2010年3月23日 (火)

確率と分布関数(10)(線形回帰の基礎)

確率と分布関数の続きです。 

 まず,最小二乗法の話をします。これは必ずしも確率の話ではなく,実際の観測値をモデル予測式を想定して予測する手法の1つです。

すなわち,多数の実測値データとモデル式による予測値の差の二乗和を最小にする予測モデル式中の最適パラメータ値を求める方法です。

[定義13-1]:x-y平面上にn個の点(xj,yj)(j=1,2,..,n)(散布図)が与えられているとする。

幾つかのパラメータを持つモデル式:y=f(x)に対し,x=xjにおける実測値yjと予測値y=f(xj)の差εj≡yj-f(xj)の二乗和:Q≡ε12+ε22+..+εn2=Σj=1n[yj-f(xj)]2を最小にする曲線(=回帰曲線:regression curve)f(x)(またはそのパラメータ)を定める方法を最小二乗法(method of minimum square)という。

 特に,予測式:f(x)が直線f(x)=β1+β2xの場合はx=xjにおける差はεj≡yj-(β1+β2j)(j=1,2,..,n)である。二乗和:Q=Q(β12)≡ε12+ε22+..+εn2=Σj=1n[yj-(β1+β2j)]2を最小にするβ12の値を求める方法は特に線形最小二乗法(linear method)といわれる。

線形最小二乗法でQの最小値を与える係数値をβ1=β^12=β^2とするとき,最適直線y=β^1+β^2xを回帰直線(regression line),係数β^1,β^2を回帰係数(regression coefficients)と呼ぶ。

12)が最小になるための2つの必要条件は∂Q/∂β1=-2Σj=1n[yj-(β1+β2j)]=0 ,および∂Q/∂β2=-2Σj=1n[xj{yj-(β1+β2j)}]=0 である。これを正規方程式(normal equations)という。

正規方程式を満たす解の値β1=β^12=β^2を解けば,β^2={Σj=1njj-n<x><y>}/{Σj=1nj2-n<x>2}=Σj=1n(xj-<x>)(yj-<y>)/{Σj=1n(xj-<x>)2},β^1=<y>-β^2<x>を得る。

 

ただし,<>は標本平均,すなわち<x>≡Σj=1nj/n,<y>≡Σj=1nj/nである。

さらに,∂Q2/∂β12=2n,∂Q2/∂β1∂β2=2Σj=1nj,∂Q2/∂β22=2Σj=1nj2なので,(∂Q2/∂β1∂β2)2-(∂Q2/∂β12)(∂Q2/∂β22)=4(Σj=1nj)2-4n(Σj=1nj2)=-4nΣj=1n (xj2-<x>)2≦0 が成り立つ。

 

それ故,β1=β^12=β^2は確かにQの"唯一の極小値=最小値"を与える。

[定理13-2]:ある変量Yが他の変数xに対してY=β1+β2x+ε (εは正規分布:N[0,σ2]に従う誤差)と表わされるとき,最尤法におけるβ12の推定値は最小二乗法による回帰係数と一致する。

(証明):x=xjのときのYの測定値をyjとすると,Yは正規分布N[β1+β2j2]に従うため,"Y1,Y2,..,Ynの同時確率分布関数=尤度関数"はL=(2πσ2)-n/2exp[-Σj=1n{yj-(β1+β2j)}2/(2σ2)]=(2πσ2)-n/2exp[-Q(β12)/(2σ2)]で与えられます。

 したがって,尤度(likelihood)Lを最大にするβ12がQ=Q(β12)を最小にするβ12に等しいことは明らかです。(証明終わり)

[定理13-3]:[定理13-2]と同じ条件下でY=β1+β2x+εを与える推定量β1=β^12=β^2は不偏推定量(unbiased estimator)である。

(証明):今のケースではβ^2=Σj=1n(xj-<x>)(Yj-<Y>)/{Σj=1n(xj-<x>)2},β^1=<Y>-β^2<x>です。

1,x2,..,xnが確定値であるのに対し,これらに対応するY1,Y2,..,Ynは様々な値を取る母集団の任意標本であると考えられます。

すると,Yj=β1+β2j+εj,かつE[εj]=0 からE[β^2]=Σj=1n(xj-<x>)E[Yj-<Y>]/{Σj=1n(xj-<x>)2}=Σj=1n(xj-<x>){β2(xj-<x>)}/{Σj=1n(xj-<x>)2}=β2,E[β^1]=E[<Y>-β^2<x>]=β1+β2<x>-β2<x>=β1を得ます。(証明終わり)

[定理13-4]:回帰係数β^1,β^2の分散,共分散は次の性質を有する。

 

(1)Var[β^1]=σ2j=1nj2)/{nΣj=1n(xj-<x>)2},Var[β^2]=σ2/{Σj=1n(xj-<x>)2},Cov[β^1,β^2]=-σ2<x>/{Σj=1n(xj-<x>)2}である。

 

(2)V≡Σj=1n{Yj-(β1+β2j)}2/(n-2)は,"分散:σ2の不偏推定量=不偏分散"である。

(証明):(1)β^2=Σj=1n(xj-<x>)(Yj-<Y>)/{Σj=1n(xj-<x>)2}にYj=β1+β2j+εj,および<Y>=β1+β2<x>+<ε>を代入すると,β^2=Σj=1n(xj-<x>){β2(xj-<x>)+εj-<ε>}/{Σj=1n(xj-<x>)2}となります。

 そこで,β^2-β2=Σj=1n(xj-<x>)(εj-<ε>)/{Σj=1n(xj-<x>)2}ですが,E[β^2-β2]=0 なのでVar[β^2]=E[(β^2-E[β2])2]=E[(β^2-β2)2]=E[{Σj=1n(xj-<x>)(εj-<ε>)}2]/{Σj=1n(xj-<x>)2}2を得ます。

ところが,E[εj]=E[<ε>]=0,E[εj]=σ2なのでVar[β^2]=E[Σj=1n(xj-<x>)2j-<ε>)2]/{Σj=1n(xj-<x>)2}2=Σj=1n(xj-<x>)2E[(εj-<ε>)2]/{Σj=1n(xj-<x>)2}2=E[εj2]/{Σj=1n(xj-<x>)2}=σ2/{Σj=1n(xj-<x>)2}です。

それ故,Var[β^1]=Var[Y-β^2<x>]=Var[Y]+<x>2Var[β^2]=σ2/n+σ2<x>2/{Σj=1n(xj-<x>)2}=σ2j=1nj2)/{nΣj=1n(xj-<x>)2}が得られます。

また,Cov[β^1,β^2]=E[β^1β^2]-β1β2=E[(<Y>-β^2<x>)β^2]-β1β2=E[<Y>β^2]-(Var[β^2]+β22)<x>-β1β2です。

ここで<Y>=β1+β2<x>+<ε>,E[β^2]=β2より,E[<Y>β^2]=β1β2+β22<x>+E[β^2<ε>]です。

右辺最後の項は,E[β^2<ε>]=E[Σj=1n(xj-<x>)(Yj-<Y>)<ε>]/{Σj=1n(xj-<x>)2}です。

 

しかし,Yj-<Y>=β2(xj-<x>)+εj-<ε>より,右辺の分子=E[Σj=1n(xj-<x>)(Yj-<Y>)<ε>]=Σj=1n(xj-<x>)E[(Yj-<Y>)<ε>]=β2E[<ε>]Σj=1n(xj-<x>)2+Σj=1n(xj-<x>)E[(εj-<ε>)<ε>]=0 です。

故に,結局Cov[β^1,β^2]=β1β2+β22<x>-(Var[β^2]+β22)<x>-β1β2=-<x>Var[β^2]=-σ2<x>/{Σj=1n(xj-<x>)2}を得ます。

(2)<Y>=β^1+β^2<x>+<ε>なので,Yj-β^1-β^2j=Yj-<Y>-β^2(xj-<x>)です。

そこでj=1n(Yj-β^1-β^2j)2=Σj=1n(Yj-<Y>)2+Σj=1nβ^22(xj-<x>)2-2β^2Σj=1n(Yj-<Y>)(xj-<x>)=Σj=1n(Yj-<Y>)2+Σj=1nβ^22(xj-<x>)2-2β^2Σj=1n(xj-<x>){β^2(xj-<x>)+εj-<ε>}です。

 

つまり,Σj=1n(Yj-β^1-β^2j)2=Σj=1n(Yj-<Y>)2-Σj=1nβ^22(xj-<x>)2-2β^2Σj=1n(xj-<x>)(εj-<ε>)と書けます。

故に,E[Σj=1n(Yj-β^1-β^2j)2]=Σj=1nE[(Yj-<Y>)2]-Σj=1n(xj-<x>)2E[β^22]-2Σj=1n(xj-<x>)E[β^2j-<ε>)]=Σj=1nE[Yj2]-nE[<Y>2]-Σj=1n(xj-<x>)2E[β^22]なる式を得ます。

ところが,E[Yj2]=Var[Yj]+E[Yj]2=σ2+(β1+β2j)2,E[<Y>2]=Var[<Y>]+E[<Y>]2=σ2/n+(β1+β2<x>)2j=1n(xj-<x>)2E[β^22]=Σj=1n(xj-<x>)2(Var[β^2]+E[β^2]2)=σ2+β22Σj=1n(xj-<x>)2です。

 したがって,E[Σj=1n(Yj-β^1-β^2j)2]=nσ2+Σj=1n1+β2j)2-σ2-n(β1+β2<x>)2-σ2+β22Σj=1n(xj-<x>)2=(n-2)σ2となります。

以上から,E[Σj=1n{Yj-(β^1+β^2j)}2/(n-2)]=σ2が示されました。(証明終わり)

さて,これまでの論議は線形最小二乗法において,変量Y=β1+β2x+εなる式の独立変数がx個だけのいわゆる単回帰です。

以下では,これをY=x1β1+x2β2+..+xpβp+ε=Σk=1pkβk+εなのように独立変数がp個の多変数x1,x2..,xpの重回帰モデルに拡張します。(ただし,x1β1が単に定数β1があるようなx=1(一定)の式も含みます。)

[定義13-5]:n個の標本Y1,Y2,..,YnをYj=xj1β1+xj2β2+..+xjpβp+εj=Σk=1pjkβk+εj (j=1,2,..,n)と置く。ただしβk(k=1,2,..,p)は未知母数であり,xjk(j=1,2,..,n;k=1,2,..,p)は既知定数である。

εj(j=1,2,..,n)は誤差で,これらは互いに独立な確率変数であって全てN[0,σ2]に従うと仮定する。これを線形回帰モデルという。

β^kをβkの推定量とするとき,Yj-Σk=1pjkβ^k (j=1,2,..,n)を残差(residual)という。残差の平方和:Q(β^)=Σj=1n(Yj-Σk=1pjkβk^)2を最小にする回帰係数:β^=(β^1,β^2,..,β^p)を求める方法を最小二乗法という。

先にQ(β^1,β^2)の最小値を与えるβ^1,β^2を係数とする回帰直線y=β^1+β^2xを求めた手続きを単回帰と呼ぶのに対して,上記のp≧3の手続きを重回帰(multiple linear regression)という。

そして,Q(β^)が最小になるための必要条件:-(1/2)(∂Q/∂β^k)=Σj=1n[xjk{Yj-(xj1β^1+xj2β^2+..+xjpβ^p)}]=0 (k=1,2,..,p)を正規方程式と呼ぶ。

 

※(注):具体的に正規方程式を解きます。

 

正規方程式はΣj=1n(xjkj1β^1+xjkj2β^2+..+xjkjpβ^p)=Σj=1njkj (k=1,2,..,p)です。

kl≡Σj=1njkjl,Dk≡Σj=1njkjと置けば,正規方程式はΣl=1pklβ^l=Dk (k=1,2,..,p)と書けます。これはp×pの係数行列Sとp次元縦ベクトル:β^≡t(β^1,β^2,..,β^p),t(D1,D2,..,Dp)による行列表現ではSβ^=です。

また,t(Y1,Y2,..,Yn),kt(x1k,x2k,..,xnk) (k=1,2,..,p)なるn次元縦ベクトルを用いるとSkl=Σj=1njkjltkl,Dk=Σj=1njkjtk(k,l=1,2,..,p)ですから,n×p行列XをX≡(1,2,..,p)で定義すると,S=tXX,tと書けます。

さらに,Yj=Σk=1pjkβk+εj (j=1,2,..,n)もβt12,..,βp),εt12,..,εn)により,=Xβεですからt=Sβtεです。

S=tXXについてdetS≠0 と仮定すれば,T=(tij)≡S-1が存在するので正規方程式:Sβ^=から,"β^=S-1=Tなる解=回帰係数"が得られます。

X=(1,2,..,p)は既知定数成分の行列なのでS=tXX,T=S-1も確率変数ではないため,E[β^]=S-1E[]=S-1E[Sβtε]=E[β]+S-1 tXE[ε]=βよりβ^はβ^の不偏推定量です。(注終わり)※

[定理13-6]:Cov(β^i,β^j)=tijσ2,特にVar(β^i)=tiiσ2 (i,j=1,2,..,p)である。またQ(β^)=Σj=1nj2-Σk=1pkβ^kであり,E[Q]=(n-p)σ2である。

(証明):Cov[Di,Dj]=Cov[Σk=1nkikl=1nljl],Cov[Yk,Yl]=Cov[εkl]=δklσ2,故にCov[Di,Dj]=Σk=1nkikjσ2=Sijσ2です。

これとβ^=S-1=TによってCov[β^i,β^j]=(S-1)ik(S-1)jlklσ2=(S-1)ijσ2を得ます。すなわち,Cov[β^i,β^j]=tijσ2,特にj=iならVar[β^i]=tiiσ2 (i,j=1,2,..,p)です。

また,Q(β^)=Σj=1n(Yj-Σk=1pjkβ^k)2=Σj=1nj2-2Σj=1nΣk=1pjkjβ^k+Σj=1nk=1pjkβ^k)2=Σj=1nj2-2Σk=1pkβ^k+Σk,l=1pklβ^kβ^l=Σj=1nj2-2Σk=1pkβ^k+Σk=1pkβ^k=Σj=1nj2-Σk=1pkβ^kと書けます。

また,Uj≡Yj-Σk=1pjkβ^kと置けばQ(β^)=Σj=1nj2であり,σ2=Var[Yj]=Var[Σk=1pjkβ^k]+Var[Uj]=Σk,l=1pjkjlCov[β^k,β^l]+Var[Uj]です。

 

したがって,Var[Uj]=σ2-Σk,l=1pjkjlklσ2を得ます。そしてE[Uj]=E[Yj-Σk=1pjkβ^k]=0です。

それ故,E[Q]=E[Σj=1nj2]=Σj=1nE[Uj2]=Σj=1n Var[Uj]=nσ2-Σj=1nΣk,l=1pjkjlklσ2=nσ2-Σk,l=1pklklσ2=nσ2-Σk=1pσ2=(n-p)σ2です。(証明終わり)

今日はここまでにします。(つづく)

参考文献:藤沢武久 著「新編 確率・統計」(日本理工出版会)

     

 

ブックオフオンライン 

iconオンライン書店 boople.com(ブープル)

|

« 最近考えていること(場の理論等覚え書き) | トップページ | 確率と分布関数(11)(区間推定)(終了) »

309. 確率・統計」カテゴリの記事

コメント

コメントを書く



(ウェブ上には掲載しません)




トラックバック


この記事へのトラックバック一覧です: 確率と分布関数(10)(線形回帰の基礎):

« 最近考えていること(場の理論等覚え書き) | トップページ | 確率と分布関数(11)(区間推定)(終了) »