統計科学への応用 統計科学への応用

1.

回帰モデル（誤差項が正規分布やｔ分布）

2.

トービットモデル（打ち切りのある回帰モデル）

3.

プロビットモデル

4.

ロジットモデル

5.

見かけ上無関係な回帰モデル

1. 1. 回帰モデル 回帰モデル

線形回帰モデル

• y _i ：確率変数で被説明変数 , i=1,…,n

• x _i =(x _i1 ,…,x _iK ) ´ : 既知の定数で説明変数

• β =( β ₁ ,…, β _K )‘: 未知のパラメータ

• ε _i 〜 N(0, σ ² ) は確率変数で誤差項 .

1. 1. 回帰モデル（続き） 回帰モデル（続き）

• 回帰モデルは

と表現でき、尤度関数は

• 事前分布

1. 1. 回帰モデル（続き） 回帰モデル（続き）

•

事後分布

ただし

1. 1. 回帰モデル（続き） 回帰モデル（続き）

（例）気温

(Y)

とビール・発泡酒の課税数量（千キロ リットル）

(X)

の関係

300 500 700 900

0 10 20 30

ビー ル

・ 発 泡 酒

( 千

ｋ ｌ︶

1. 1. 回帰モデル（続き） 回帰モデル（続き）

例： WinBUGS のコード

•

モデルの指定

model{

# ｆ（ｘ|β,σ²) τ=σ-² ← 尤度関数

for(i in 1:n){ #nは観測値の個数

mu[i] <-beta[1]+beta[2]*x[i] # μ_i=β₁+β₂x_i y[i]~dnorm(mu[i],tau) # y_i ~N(μ_i,σ²) }

# π(β,σ²) ← 事前分布

for(i in 1:K){ #Kはベータの次元=2 beta[i]~dnorm(0,1.0E-6)} #β_i~N(0,10⁶)

tau~dgamma(1.0E-6,1.0E-6)

1. 1. 回帰モデル（続き） 回帰モデル（続き）

•

計算結果のまとめ

584.7 429.9

271.3 0.833

78.34 429.2

β

₁

28860 10150

4593 78.13

6550 11760

σ

²

17.42 8.763

0.286 0.045

4.305 8.813

β

₂

97.5 %

中央

値

2.5

％

MC error

標準 偏差 事後

平均

1. 1. 回帰モデル（続き） 回帰モデル（続き）

beta[1] sample: 10000

-250.0 0.0 250.0 500.0 750.0

beta[2] sample: 10000

-20.0 0.0 20.0 40.0

sigma2 sample: 10000

1. 1. 回帰モデル（続き） 回帰モデル（続き）

• y _n+1 の予測分布は・・・

• ギブス・サンプラーの最後のステップで , 得ら れた ( β , σ ² ) のもとで y _n+1 | β , σ ² , x _n+1

~N(x _n+1 β , σ ² ) を発生させればよい .

1. 1. 回帰モデル（続き） 回帰モデル（続き）

誤差項がｔ分布であるような回帰モデル

–

標準的な回帰モデルでは

,

誤差項は ε

_i ~N(0,

σ

² )

と仮定した

.

–

これを正規分布ではなく

, t

分布にしたモデルも 容易に扱うことができる

.

とおくと（ν

₀

は既知の定数）、ε

_i

の周辺分布は自

1. 1. 回帰モデル（続き） 回帰モデル（続き）

• 回帰モデルは

と表現でき、尤度関数は

• 事前分布

1. 1. 回帰モデル（続き） 回帰モデル（続き）

•

事後分布

ただし

2. 2. トービットモデル トービットモデル

トービットモデル（打ち切りのある回帰モデル）

•

被説明変数

y

_iが

,

ある閾値を超えたときのみに観測される 回帰モデル

.

閾値を

0

とすれば

,

• y

_i^*：確率変数で被説明変数

, i=1,…,n

• y

_i：

y

_i^*

>0

のとき観測され

,

そうでないとき観測されない

.

• x

_i

=(x

_i1

,…,x

_iK

)

´

:

既知の定数で説明変数

•

β

=(

β₁

,…,

β_K

)‘:

未知のパラメータ

•

ε_i 〜

N(0,

σ²

)

は確率変数で誤差項

.

2. 2. トービットモデル（続き） トービットモデル（続き）

• トービットモデルは

と表現でき、尤度関数は

• 事前分布

2. 2. トービットモデル（続き） トービットモデル（続き）

•

事後分布

ただし

2. 2. トービットモデル（続き） トービットモデル（続き）

例：既婚女性の 1 年間の労働時間 (y) を

1.

子供の人数

(x1), 2.

年齢

(x2),

3.

教育年数

(x3), 4.

夫の年収

(x4)

で説明するモデルを考える

.

–

労働時間は、就職していない人は

0

時間になる ので、トービットモデルを用いる。

2. 2. トービットモデル（続き） トービットモデル（続き）

計算例（Ox言語）

for(cn=-cburn;cn<crepeat;++cn){

mB1=invert((1/dsigma2)*mx‘*mx+invert(mB0)); #βのサンプリング vb1=mB1*((1/dsigma2)*mx'*vystar+invert(mB0)*vb0);

vb=vb1+choleski(mB1)*rann(cm,1);

for(i=0;i<cnobs;++i){ #y*のサンプリング if(vy[i]==0){

do{vystar[i]=mx[i][]*vb+rann(1,1)*sqrt(dsigma2);

}while(vystar[i] >=0);

}

else{vystar[i]=vy[i];}}

＃σ²のサンプリング

da=0.5*(0.002+cnobs);db=0.5*(0.002+(vystar-mx*vb)'*(vystar-mx*vb));

dsigma2=1/rangamma(1,1,da,db);

}

2. 2. トービットモデル（続き） トービットモデル（続き）

計算結果

49.01 -24.03

18.77 12.49

β

₁

0.56 -13.38

3.53 -6.15

β

₂

-0.30 -1.66

0.35 -0.95

β

₃

6.38 0.70

1.45

β

₄ 3.46

5.44E-4 -2.31E-4

1.95E-4 1.57E-4

β

₅

97.5%

2.5

％ 標準偏差

事後平均

2. 2. トービットモデル（続き） トービットモデル（続き）

3. 3. プロビットモデル プロビットモデル

プロビットモデルー選択行動のモデル

• 被説明変数 y _i が 0 または 1 といった離散的な 値をとる .

–

例

:

新車を購入する

(y _i =1),

しない

(y _i =0)

で ある

.

–

例：ある顧客にローンを貸与する

(y _i =1),

しな い

(y _i =0)

3. 3. プロビットモデル（続き） プロビットモデル（続き）

•

モデル

–

ここでF(・)は0から1の間をとる関数

– F

に標準正規分布の分布関数→プロビットモデル

– F

にロジスティック分布の分布関数→ロジットモデル

•

潜在変数

y

_i

*

を

と定義すると

y*

が与えられたとき

,

標準的な回帰モデルで 分散が

1

と既知のモデルとなる

.

3. 3. プロビットモデル（続き） プロビットモデル（続き）

• プロビットモデルは

と表現でき、尤度関数は

• 事前分布

3. 3. プロビットモデル（続き） プロビットモデル（続き）

•

事後分布

ただし

3. 3. プロビットモデル（続き） プロビットモデル（続き）

例

:

地方学校税に関する住民投票の行動分析のために

,

有権 者に対する調査

.

そのなかの

95

人のデータが

Pindyck and Rubinfeld (1998)

にあり

,

このデータに対してプロビットモ デルをあてはめる

.

– Y:

学校税率引き上げに賛成のとき

1,

反対のとき

0 ¥¥

– X_1:

定数項

– X_2:子供が私立学校に在籍しているとき 1, いないとき 0 – X_3:

（私立または公立）学校の教師であるとき

1,

ないとき

0 – X_4:家計所得(年額）の対数値

– X_5:

財産税

(

年額

)

の対数値

3. 3. プロビットモデル（続き） プロビットモデル（続き）

例： Ox 言語の例

for(cn=-cburn;cn<crepeat;++cn){

#β~N(b₁,B₁)

vb1=mB1*(mx'*vystar+invert(mB0)*vb0);

vb=vb1+choleski(mB1)*rann(cm,1);

# y*のサンプリング

for(i=0;i<cnobs;++i){

if(vy[i]==0){do{vystar[i]=mx[i][]*vb+rann(1,1);

}while(vystar[i] >=0);}

else{do{vystar[i]=mx[i][]*vb+rann(1,1);

}while(vystar[i] <0);}

} }

3. 3. プロビットモデル（続き） プロビットモデル（続き）

3. 3. プロビットモデル（続き） プロビットモデル（続き）

3. 3. プロビットモデル（続き） プロビットモデル（続き）

計算結果

2.43 -12.3

3.73 -4.91

β

₁

0.48 -1.23

0.44 -0.39

β

₂

4.00 0.55

0.89

β

₃ 2.05

2.36 0.62

0.45

β

₄ 1.44

-0.18 -2.55

0.60 -1.33

β

₅

97.5%

2.5

％ 標準偏差

事後平均

4. 4. ロジットモデル ロジットモデル

•

ロジットモデルでは

•

事前分布

•

尤度関数

4. 4. ロジットモデル ロジットモデル ( ( 続き 続き ) )

•

事後分布

• MH

アルゴリズム

–

提案分布は正規分布

N(

ｍ

,V).

ただし

,

4. 4. ロジットモデル ロジットモデル ( ( 続き 続き ) )

• 提案分布（正規分布 , 密度関数は

φ ( β |m,V) とおく）からの候補β´を採択す る確率は ,

となる .

–

※比例定数の項は分子分母で共通するためキャ ンセルアウトする

.

4. 4. ロジットモデル ロジットモデル ( ( 続き 続き ) )

例：

WinBUGS

のコード

•

モデルの指定

model {

for (i in 1:n) { y[i] ~ dbern(p[i])

logit(p[i]) <- b[1] + b[2]*x[i,1]+ b[3]*x[i,2] + b[4]*x[i,3]

+ b[5]*x[i,4]

}

#prior

for (j in 1:K) {b[j] ~ dnorm(0,0.001)}

}

3. 3. プロビットモデル（続き） プロビットモデル（続き）

例：

Ox

言語の例

MH

の採択率は

72.75%

for(cn=-cburn;cn<crepeat;++cn){

vb_o=vb;fLK(vb_o,&dl_o,0,0); #現在の点での事後密度 dh_o=dlhat+vl1hat‘*(vb_o-vbhat) #現在の点での提案密度

+0.5*(vb_o-vbhat)'*ml2hat*(vb_o-vbhat);

vb_n=vb1+mB1root*rann(cm,1); #候補点の提案

fLK(vb_n,&dl_n,0,0); #候補点での事後密度 dh_n=dlhat+vl1hat‘*(vb_n-vbhat) #候補点での提案密度

+0.5*(vb_n-vbhat)'*ml2hat*(vb_n-vbhat);

dfrac=exp(dl_n+dh_o-dl_o-dh_n); #採択確率の計算 if(ranu(1,1)<=dfrac){vb=vb_n;}

}

4. 4. ロジットモデル ロジットモデル ( ( 続き 続き ) )

4. 4. ロジットモデル（続き） ロジットモデル（続き）

計算結果

3.87 -23.2

6.91 -9.66

β

₁

0.84 -2.12

0.76 -0.64

β

₂

6.95 0.92

1.56

β

₃ 3.55

4.18 1.00

0.81

β

₄ 2.51

-0.20 -4.30

1.06 -2.15

β

₅

97.5%

2.5

％ 標準偏差

事後平均

5. 5. 見かけ上無関係な回帰モデル 見かけ上無関係な回帰モデル

Seemingly Unrelated Model

• y _i

：

mx1

確率変数ベクトルで被説明変数

, i=1,…,n

–

例：

t

時点における

m

個の危険資産の収益率など

• X _i =

既知の

mxK

で説明変数行列

•

β

=(

β

₁ ,…,

β

_K )‘:

未知のパラメータ

•

ε

_i

〜

N(0,

Σ

)

はｍｘ１確率変数ベクトルで誤差項

.

5.SUR

5.SUR モデル（続き） モデル（続き）

• 事前分布

• 尤度関数は

5.SUR

5.SUR モデル（続き） モデル（続き）

• 事後分布

• ただし

5.SUR

5.SUR モデル（続き） モデル（続き）

例 :General Electric と Westinghouse

– Gross investment (Y):

被説明変数

– Value of the Firm (V)

：説明変数

– Stock of plant and equipment(K) – y _i =

β

₁ +

β

₂ V _i

＋β

₃ K _i

＋ε

_i

– Grunfeld Investment Data, 100 yearly

observations. 1935-1954.

5.SUR

5.SUR モデル モデル ( ( 続き 続き ) )

例：

WinBUGS

のコード

•

モデルの指定

model {

for (t in 1:T) {

y[t,1:2] ~ dmnorm(mu[t,],Omega[ , ]);

mu[t,1] <- beta[1,1]+beta[1,2]*v[t,1]+beta[1,3]*k[t,1]

mu[t,2] <- beta[2,1]+beta[2,2]*v[t,2]+beta[2,3]*k[t,2]

}

# priors on regression coefficients

for (i in 1:M) {for (j in 1:3) {beta[i,j] ~ dnorm(0,0.001)}}

Omega[1 : M , 1 : M] ~ dwish(R[ , ], 2)

# cross-correlation matrix of dimension M=2

5.SUR

5.SUR モデル モデル ( ( 続き 続き ) )

beta[1,1] sample: 50000

-200.0 -100.0 0.0 100.0

beta[1,2] sample: 50000

-0.05 0.0 0.05 0.1

beta[1,3] sample: 50000

-0.1 0.0 0.1 0.2 0.3

rho sample: 50000

-0.5 0.0 0.5 1.0 1.5

5.SUR

5.SUR モデル モデル ( ( 続き 続き ) )

beta[2,1] sample: 50000

-50.0 -25.0 0.0 25.0 50.0

beta[2,2] sample: 50000

-0.05 0.0 0.05 0.1 0.15 beta[2,3] sample: 50000

5. 5. SUR SUR モデル（続き） モデル（続き）

•

計算結果のまとめ

0.18 0.14

0.08 4.9E-4

0.03

β

₁₃ 0.13

0.05 0.0.03

0.01 3.7E-4

0.01

β

₁₂ 0.03

14.33 1.03

-11.71 0.21

6.67

β

₂₁ 1.13

0.08 0.05

0.03 4.8E-4

0.01

β

₂₂ 0.05

26.65 -14.05

-54.88 0.65

20.78 -14.21

β

₁₁

0.90 0.76

0.46 0.001

0.12 Ρ 0.74

0.18 0.07

-0.05 0.001

0.06

β

₂₃ 0.06

97.5%

2.5

％ 中央値

MC rror

標準偏差 事後平均

ドキュメント内 untitled (ページ 57-100)