1 単回帰モデル

最尤法
宮﨑憲治
2016 年 1 月 18 日
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
1 / 32
.
目次
1
最尤法
2
ベルヌーイ分布
3
正規分布
4
正規線形回帰モデル
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
2 / 32
.
最尤法
n 個の確率変数 {xi }ni=1 の確率密度関数があるパラメータ θ に依
存するとき fn (x1 , . . . , xn ; θ) と記述する.
確率密度関数 fn (x1 , . . . , xn ; θ) を母数 θ の関数とみなすとき, 尤度
関数 L(θ) といい, 尤度を最大にするように θ を選ぶ方法を最尤法
(Method of Maximum Likelihood) という.
尤度関数を最大にする θの値と, 対数尤度関数
ℓ(θ) = log L(θ) = log fn (x1 , . . . , xn ; θ)
を最大値にする θ の値は同じなので, 数学的扱いやすい後者を用
いる.
独立同一分布にしたがっているのなら, fn = f (x1 ) × · · · × f (xn )
であり,
∑
ℓ(θ) =
log f (xi ; θ)
である.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
3 / 32
.
最尤法
最尤推定量は尤度関数を最大とする推定量である. 正則条件のも
と, 最尤推定量 θ̂ は一致で, 漸近正規である.
√
d
n(θ̂ − θ) → N (0, I(θ̂)−1 )
ここで, I(θ) はフィッシャー情報量とよばれ, 対数尤度関数 ℓ(θ) の
一次導関数 ℓ′ (θ) をもちいて,
I(θ̂) = E[{ℓ′ (θ̂)}2 ]/n
と定義される.
また正則条件のもと, 対数尤度関数の二次導関数 ℓ′′ (θ) をもちいて,
I(θ̂) = −E[ℓ′′ (θ̂)]/n
となることが知られている.
パラメータ θ 複数の場合にも拡張可能であり, その際はフィッシ
ャー情報量 I(θ̂) は行列になり, I(θ̂)−1 は逆行列である.
また漸近正規の推定量の中で分散が一番小さい, つまり漸
近有効であることが知られている.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
4 / 32
.
ベルヌーイ分布
一枚のコインを 5 回投げて「表・表・裏・表・裏」とでたとき, 表
を 1 とし, 裏を 0 とし, コインの表の出る確率 p を最尤推定法で推
定する.
ある p のもと, この事象がおこる確率を L(p) とすると以下になる:
L(p) = p × p × (1 − p) × p × (1 − p)
= p3 (1 − p)2
よって尤度関数を最大にする p が最尤法であり, p = 0.6 となる.
また, 対数尤度関数は以下になる:
ℓ(p) = 3 log p + 2 log(1 − p)
この場合でも最大値は p = 0.6 である.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
5 / 32
.
ベルヌーイ分布
curve(3*log(x)
+2*log(1-x),0,1)
-6
-8
-12
-10
3 * log(x) + 2 * log(1 - x)
0.020
0.015
0.010
0.005
0.0
0.2
0.4
0.6
0.8
-14
0.000
x^3 * (1 - x)^2
0.025
-4
0.030
0.035
curve(x^3*(1-x)^2,0,1)
1.0
0.0
x
0.2
0.4
0.6
0.8
1.0
x
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
6 / 32
.
R の実践
Be(0.5) にしたがう乱数を 100 個発生させ, 対数尤度関数を作成する.
rm(list=ls(all=TRUE))
x <- rbinom(100,1,0.5)
p<-seq(0.01,0.99,0.01)
lh<-numeric(length(p))
for (i in 1:length(p)) {
lh[i]<-sum(x*log(p[i])+(1-x)*log(1-p[i]))
}
plot(p,lh,type='l')
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
7 / 32
.
R の実践
最尤推定量を求める.
-150
## [1] -68.99438
-200
p[which.max(lh)]
## [1] 0.46
round(mean(x),2)
-250
lh
-100
max(lh)
0.0
0.2
0.4
0.6
0.8
1.0
p
## [1] 0.46
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
8 / 32
.
ベルヌーイ分布
密度関数:
f (x; p) = px (1 − p)1−x
尤度関数:
L(p) =
n
∏
f (xi ; p)
i
対数尤度関数:
log L(p) = log p
∑
xi + log(1 − p)
.
宮﨑憲治
最尤法
∑
.
.
.
.
.
(1 − xi )
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
9 / 32
.
ベルヌーイ分布
最大化のための 1 階条件:
∑
∂ℓ(p)
1∑
1
xi ) = 0
=
(n −
xi +
∂p
p
1−p
これを整理すると:
∑
p̂ =
xi
n
つまり, 最尤法は標本平均と一致する.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
10 / 32
.
ベルヌーイ分布
2 階微分:
∂ 2 ℓ(p)
∂p2
∑
1 ∑
1
(n −
xi )
xi +
2
2
p
(1 − p)
∑
(1 − 2p) xi + p2 n
p2 (1 − p)2
= −
=
その期待値:
E[
∂ 2 ℓ(p)
(1 − 2p)np + p2 n
n
]
=
=
2
2
2
∂p
p (1 − p)
p(1 − p)
フィッシャー情報量: I(p) = p(1 − p)
最尤推定量:
√
d
n(p̂ − p) → N (0, p(1 − p))
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
11 / 32
.
正規分布
平均身長を調べるため標本として 3 人抽出したところ「162, 171,
183」であった. 母集団が正規分布 N (µ, σ 2 ) にしたがうと仮定し
て, 平均身長を最尤推定法で推定する.
密度関数:
]
[
(x−µ)2
1
1
(x − µ)2
f (x) = √
e− 2σ2 = √
exp −
2σ 2
2πσ 2
2πσ 2
ある µ のもと, この事象がおこる確率を L(µ) とすると以下になる:
(162−µ)2
(171−µ)2
(183−µ)2
1
1
1
√
e− 2σ2 × √
e− 2σ2 × √
e− 2σ2
2πσ 2
2πσ 2
2πσ 2
)3
[
]
(
1
(162 − µ)2 + (171 − µ)2 + (183 − µ)2
√
=
exp −
2σ 2
2πσ 2
L(µ) =
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
12 / 32
.
正規分布
よって最尤法は
(162 − µ)2 + (171 − µ)2 + (183 − µ)2 = 3(µ − 172)2 + 222
を最大にする値で µ̂ = 172 である.
どの σ 2 をとっても µ̂ = 172 である.
対数尤度関数は
3
ℓ(µ) = log L(µ) = − (log σ 2 + log 2π)
2
(162 − µ)2 + (171 − µ)2 + (183 − µ)2
−
2σ 2
なので, 最大になる値は µ̂ = 172 で一致する.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
13 / 32
.
正規分布
N (2, 1) にしたがう乱数を 100 個発生させ, 対数尤度関数を作成する.
rm(list=ls(all=TRUE))
x <- rnorm(100,2,1)
mu<-seq(-1,3,0.01)
lh<-numeric(length(mu))
for (i in 1:length(mu)) {
lh[i]<-sum(log(dnorm(x,mu[i],1)))
}
plot(mu,lh,type='l')
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
14 / 32
.
正規分布
最尤推定量を求める.
-200
max(lh)
-400
mu[which.max(lh)]
-500
## [1] 2.02
-600
lh
-300
## [1] -157.0637
round(mean(x),2)
-1
0
1
2
3
mu
## [1] 2.02
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
15 / 32
.
正規分布
密度関数:
[
1
f (x; µ, σ ) = √
exp − 2 (xi − µ)2
2σ
2πσ 2
2
1
尤度関数:
L(µ, σ 2 ) =
n
∏
]
f (xi ; µ)
i
対数尤度関数:
n
1 ∑
ℓ(µ, σ 2 ) = log L(µ, σ 2 ) = − (log σ 2 + log 2π) − 2
(xi − µ)2
2
2σ
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
16 / 32
.
正規分布
最大化のための 1 階条件 (ω = σ 2 とする):
∂ℓ(µ, ω)
∂µ
∂ℓ(µ, ω)
∂ω
1∑
(xi − µ) = 0
ω
n
1 ∑
= −
+ 2
(xi − µ)2 = 0
2ω 2ω
=
これを整理すると:
∑
µ̂ =
xi
n
ω̂ = σ̂ 2 =
1∑
(xi − µ̂)2
n
µ の最尤推定量は標本平均と一致するが, σ 2 の最尤推定量は標本
分散でないため, 不偏推定量でないが, 一致推定量である.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
17 / 32
.
正規分布
2 階微分:
∂ 2 ℓ(µ, ω)
∂µ2
2
∂ ℓ(µ, ω)
∂µ∂ω
2
∂ ℓ(µ, ω)
∂ω 2
その期待値:
n
ω
1 ∑
(xi − µ)
ω2
n
1 ∑
−
(xi − µ)2
2ω 2 ω 3
= −
=
=
[
]
∂ 2 ℓ(µ, ω)
n
E
= − 2
2
∂µ
σ
[ 2
]
∂ ℓ(µ, ω)
E
= 0
∂µ∂ω
]
[ 2
n
n
∂ ℓ(µ, ω)
= − 2 =− 4
E
2
∂ω
2ω
2σ
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
18 / 32
.
正規分布
フィッシャー情報量:
[
2
I(µ, σ ) =
その逆行列:
2 −1
I(µ, σ )
]
1/σ 2
0
0
1/(2σ 4 )
[
=
σ2 0
0 2σ 4
]
最尤推定量:
√
d
n(µ̂ − µ) → N (0, σ 2 )
√
d
n(σ̂ 2 − σ 2 ) → N (0, 2σ 4 )
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
19 / 32
.
正規線形回帰モデル
取りうる値が実数の連続確率変数 y を考える.
x を確率変数として, 正規線形回帰モデル (Normal Linear
Regression Model) は
y = α + βx + u
とも書き表せる. y を被説明変数, x を説明変数, u を誤差項とい
う.
仮定:
(ui , xi ) は i について独立同一分布にしたがう.
u は正規分布 N (0, σ 2 ) にしたがい, x と独立である
この仮定のもと, 観測値からパラメータ α と β を推定する.
推定法として最小二乗法と最尤法が有名である.
線形回帰モデルの場合, パラメータ α と β の推定方法は最小二乗
法と最尤法も同じである.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
20 / 32
.
線形回帰モデル
(u, x) の同時密度関数 f は互いに独立なため
f (x, y) = fU (u)fX (x)
と表せる.
このときの対数尤度関数 (ω = σ 2 ):
L =
n
∏
fU (yi − α − βxi )fX (xi )
i=1
ℓ(α, β, ω) = log L =
∑
log fU (yi − α − βxi ) +
∑
log fX (xi )
推計すべきパラメータが fX に含まれていなければ, 無視しても推
計効率は変わらない.
なお最小二乗法は
∑
(yi − α − βxi )2
を最小にするように α と β を求める方法である.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
21 / 32
.
線形回帰モデル
尤度関数 (ω = σ 2 ):
[
]
1
1
2
√
L=
exp − (yi − α − βxi )
2ω
2ωπ
i=1
n
∏
対数尤度関数:
n
1 ∑
ℓ(α, β, ω) = log L = − (log ω + log 2π) −
(yi − α − βxi )2
2
2ω
1 階条件:
∂ℓ(α, β, ω)
∂α
∂ℓ(α, β, ω)
∂β
∂ℓ(α, β, ω)
∂ω
1∑
(yi − α − βxi ) = 0
ω
1∑
=
(yi − α − βxi )xi = 0
ω
n
1 ∑
= −
+ 2
(yi − α − βxi )2 = 0
2ω 2ω
=
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
22 / 32
.
線形回帰モデル
最尤推定量:
∑
β̂ =
(xi − x̄)(yi − ȳ)
∑
(xi − x̄)2
α̂ = ȳ − β̂ x̄
1∑
ω̂ = σ̂ 2 =
(yi − α̂ − β̂xi )2
n
p
p
p
最尤推定量は一致推定量: α̂ → α, β̂ → β, σ̂ 2 → σ 2 .
最尤推定量は不偏推定量: E[α̂] = α, E[β̂] = β
E[σ̂ 2 ] = σ 2 (n − 2)/n ̸= σ 2 なので, 誤差項の分散の不偏推定量
SER =
∑
(yi − α̂ − β̂xi )2
n−2
を採用する.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
23 / 32
.
線形回帰モデル
最尤推定量の (x が与えられたもとの条件付き) 分布
∑
√
d
n(α̂ − α) → N [ 0, σ 2 (1 + nx̄2 /
(xi − x̄)2 )]
∑
√
d
n(β̂ − β) → N [ 0, nσ 2 /
(xi − x̄)2 ]
√
d
n(σ̂ 2 − σ 2 ) → N [ 0, 2σ 4 ]
正確にはティー分布などにしたがうが, n を十分大きく取れば正規
分布に分布収束する.
最尤推定量の (x が与えられたもとの条件付き) 共分散:
∑
Cov[α̂, β̂] = −x̄/
(xi − x̄)2
Cov[α̂, σ̂ 2 ] = Cov[β̂, σ̂ 2 ] = 0
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
24 / 32
.
多変量線形回帰モデル
線形モデルの説明変数を k 個に拡張できる:
y = β0 + β1 x1 + . . . + βk xk + u
u ∼ N (0, σ 2 )
u は正規分布 N (0, σ 2 ) にしたがい, (x1 , . . . , xk ) と独立とする.
このときも単回帰のときと同様に係数についての最尤推定量と最
小自乗推定量は同じであり, 不偏で一致推定量で漸近分散である.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
25 / 32
.
モデル比較
線形モデルの比較の際, 決定係数もしくは修正済み決定係数がよ
く使われる.
一方, 最尤推定量の場合, AIC (赤池情報量基準) が有名.
−2nℓ̂ + 2k = −2n log L̂ + 2k
n: 観測値数, k: 説明変数の数, log L̂: 最大対数尤度
より小さい値が望ましい.
ほかにも BIC (ベイズ情報量基準) というのもある.
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
26 / 32
.
R によるコード実行
set.seed(123)
n <- 300
alpha <- 3
beta <- 1
x <- runif(n)
y <- alpha + beta * x+rnorm(n,0,sd=0.1)
df <- data.frame(y = y, x = x)
上記のようにデータを生成して, 以下のコードを実行する.
fm <- lm(y ~ x, data = df)
summary(fm)
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
27 / 32
.
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
Call:
lm(formula = y ~ x, data = df)
Residuals:
Min
1Q
Median
-0.26570 -0.06375 -0.00183
3Q
0.06242
Max
0.32267
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.00872
0.01168 257.64
<2e-16 ***
x
0.98728
0.02039
48.43
<2e-16 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.09893 on 298 degrees of freedom
Multiple R-squared: 0.8873,^^IAdjusted R-squared: 0.8869
F-statistic: 2345 on 1 and 298 DF, p-value: < 2.2e-16
. . . .... .... .... . . . . .
R によるコード実行
最尤法にもとづくより一般的な回帰分析を実行するために glm という
コマンドがある.
fm <- glm(y ~ x, data = df, family=gaussian)
summary(fm)
多変量の場合も同様に次のように実行すればよい.
z <- runif(n)
dfz <- cbind(df,z)
fmz <- glm(y ~x+z, data=dfz, family=gaussian)
summary(fmz)
.
宮﨑憲治
最尤法
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
.
.
2016 年 1 月 18 日
.
.
.
.
.
.
.
29 / 32
.
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
Call:
glm(formula = y ~ x, family = gaussian, data = df)
Deviance Residuals:
Min
1Q
Median
-0.26570 -0.06375 -0.00183
3Q
0.06242
Max
0.32267
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.00872
0.01168 257.64
<2e-16 ***
x
0.98728
0.02039
48.43
<2e-16 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.009786547)
Null deviance: 25.8668
Residual deviance: 2.9164
AIC: -532.67
on 299
on 298
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 2
. . . .... .... .... . . . . .
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
##
Call:
glm(formula = y ~ x + z, family = gaussian, data = dfz)
Deviance Residuals:
Min
1Q
Median
-0.26108 -0.06284 -0.00161
3Q
0.06457
Max
0.32174
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.00324
0.01558 192.795
<2e-16 ***
x
0.98738
0.02041 48.371
<2e-16 ***
z
0.01101
0.02070
0.532
0.595
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.009810149)
Null deviance: 25.8668
Residual deviance: 2.9136
AIC: -530.95
on 299
on 297
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 2
. . . .... .... .... . . . . .
Table:
Dependent variable:
y
x
(1)
(2)
0.987∗∗∗
(0.020)
0.987∗∗∗
(0.020)
z
0.011
(0.021)
Constant
3.009∗∗∗
(0.012)
3.003∗∗∗
(0.016)
Observations
Log Likelihood
Akaike Inf. Crit.
300
268.334
−532.668
300
268.477
−530.953
Note:
∗ p<0.1; ∗∗ p<0.05; ∗∗∗ p<0.01
. . . .... .... .... . . . . .

Download Report