変な形の尤度関数をプロットする - ジョンとヨーコのイマジン日記

とりあえずプロット
R のコード
参考文献
最尤推定の基礎

とりあえずプロット

例1: 混合正規分布

モデルとして次の分布を考える.

$\displaystyle p(y; a, b) = (1-a)\mathcal{N}(y; 0,1) + a\mathcal{N}(y; 0,b)$

（この例は『ベイズ統計の理論と方法』に出てくる.）

ここで $\mathcal{N}(y; 0,1)$ は平均0, 分散1 の正規分布の密度関数である. $0 \le a \le 1$ とする.

この手のモデルは2つの分布が混ざったような形になるので混合分布とよばれ, モデルに基づくクラスタリングなどの際によく使われる.

サンプルを生成した分布（ $q(y)$ とする）が $\mathcal{N}(y,0,1)$ のとき, $a=0$ なら $b$ はなんでも, また $b=0$ なら $a$ はなんでも, データを生成した分布とモデルが一致する.

つまりサンプルを生成した分布に対して, モデルの最適なパラメータが一意に定まらない.

サンプルの大きさを100として, このときの尤度関数を等高線でプロットしてみよう.

グリッドサーチ（範囲を限定して適当な幅で単純な総当り）で求めた最尤推定値をバツ印でプロットしてあるが, これはサンプルのあらわれ方によって大きく揺らぐ.

次に $\mathcal{N}(y,0,1)$ で生成したサンプル $\mathcal{N}(y, b, 1)$ で生成したサンプルを50個ずつ混ぜてみよう

このとき, モデルの最適なパラメータは $a=0.5, b=b$ である.

$b=1$ のとき:

まずはデータのヒストグラムを見る.

まあ, 言われたらなんとなく2つ山があることに気づくかな、くらいである.

このときの尤度関数を見る.

尤度関数のモードを中心とした2変量正規分布からの乱数をグレーの点で重ねてある.

（2変量正規分布の密度関数の等高線は楕円になる 2変量正規分布の等高線（確率楕円）の描き方 - ジョンとヨーコのイマジン日記）.

モードの周りはいいとして, 裾のほうは正規分布で近似するのが苦しい感じである.

$b=0.5$ のとき:

データのヒストグラムを見る.

2つの分布が混ざっていることにはなかなか気づかないだろう.

このときの尤度関数を見る.

正規分布で近似するのは苦しい感じである.

$b=0.25$ のとき:

2つの分布が混ざっているのではないかと感じる人はほとんどいなさそうである.

このときの尤度関数を見る.

正規分布による近似は苦しい.

例2: 非線形回帰

モデルとして次の分布を考える.

$\displaystyle p(y; a, b, x) = \mathcal{N}(y, b\tanh(a x),1)$

$x$ は既知の数とする.

（この例は藤原香織、渡辺澄夫（2006）ベイズ尤度比検定による変化点検出のシミュレーション - ジョンとヨーコのイマジン日記と同じ. ）

例えば $x$ を時間として, なんらかの施策の前後でデータの変化を見たいときは（データを適当にスケーリングするとして）このようなモデルを考えるかもしれない.

ちなみに, 一番単純なニューラルネットワークではこのような非線形の変換をベクトルに対して繰り返して行い, 例えば $\mathcal{N}( y; \tanh (\tanh(x w_1) w_2) w_3, 1)$ のようなモデルを考える.

さて, $q(y) = \mathcal{N}(y; 0,1)$ としてサンプルを生成してプロットすると次のようになる.

このとき, $a=0$ のとき $b$ はなんでも, また $b=0$ のとき $a$ はなんでも, データを生成した分布とモデルが一致する.

例1のときと同様, サンプルの大きさを100として, このときの尤度関数を等高線でプロットしてみよう.

グリッドサーチで求めた最尤推定値をバツ印でプロットしてあるが, これはサンプルのあらわれ方によって大きく揺らぐ.

$x \le 0$ のとき:

$\displaystyle q(y; x) =\mathcal{N}(y; -0.4,1)$

$x > 0$ のとき:

$\mathcal{N}(y; 0.4, 1)$

で生成したサンプルを50個ずつ混ぜてみよう. （といっても $x=0$ は今回の例では登場しない.）

このとき $a = \infty, b=0.4$ または $a=-\infty, b=0.4$ でデータを生成した分布とモデルが一致する.

無限大まではプロットできないので $-5 \le a \le 5$ の $-5 \le b \le 5$ の範囲で尤度関数の等高線を描いてみる.

例3: ガンマ分布

最後に, 形状パラメータ2, 尺度パラメータ1のガンマ分布でサンプルを100個生成して, モデルを正規分布 $\mathcal{N}(y; a,b)$ とした例をやってみよう.

（この例は『In All Likelihood』に出てくる）

尤度関数のモードを中心とした2変量正規分布からの乱数をグレーの点で重ねてあるが, これは近い形の分布になっている.

（途中計算は省略して） $a=2$ , $b=2$ で, データを生成した分布とモデルの間のカルバック・ライブラ距離は最小となる.

最小となってもぴったり一致はしないので「選んだモデルの範囲内では」というただし書きつきではあるが, 例3は最尤法がうまくいっている例である.

R のコード

library(ggplot2)

#対数尤度関数
llmixnorm <- function(par,y){
  a0 <- par[1]
  b0 <- par[2]
  sum(log((1-a0)*dnorm(y)+a0*dnorm(y,b0)))
}

lltanh <- function(par,y,x){
  a0 <- par[1]
  b0 <- par[2]
  sum(dnorm(y, b0*tanh(a0*x), log=TRUE))
}

llnorm <- function(par,y, lp=FALSE){
  a0 <- par[1]
  b0 <- par[2]
  sum(dnorm(y, a0, b0, log=TRUE))
}

# 2変量正規乱数の生成
sample_norm <- function(optp, H, np){
  U <- chol(-H) #コレスキー分解
  s <- backsolve(U, matrix(rnorm(2*np),2,np))+optp # Uの逆行列を左からかけてる
  data.frame(a=s[1,], b=s[2,])
}

### Example1
N <- 100L
set.seed(1)
y0 <- rnorm(N)
y11 <- c(rnorm(N/2),rnorm(N/2,1))
y12 <- c(rnorm(N/2),rnorm(N/2,0.5))
y13 <- c(rnorm(N/2),rnorm(N/2,0.25))

ggplot()+
  geom_histogram(data=NULL, aes(x=y11), fill="grey70", bins=25)+
  theme_bw(14)+labs(x="y (1)", y="count")
ggsave("./Desktop/y1.png")

ggplot()+
  geom_histogram(data=NULL, aes(x=y12), fill="grey70", bins=25)+
  theme_bw(14)+labs(x="y (2)", y="count")
ggsave("./Desktop/y2.png")

ggplot()+
  geom_histogram(data=NULL, aes(x=y13), fill="grey70", bins=25)+
  theme_bw(14)+labs(x="y (3)", y="count")
ggsave("./Desktop/y3.png")

parms <- expand.grid(a=seq(0,1,length.out = 200),
                     b=seq(-5,5,length.out = 200))

l10 <- apply(parms, 1, llmixnorm, y=y0)
l11 <- apply(parms, 1, llmixnorm, y=y11)
l12 <- apply(parms, 1, llmixnorm, y=y12)
l13 <- apply(parms, 1, llmixnorm, y=y13)

optp10 <- parms[which.max(l10),]
optp11 <- unlist(parms[which.max(l11),])
optp12 <- unlist(parms[which.max(l12),])
optp13 <- unlist(parms[which.max(l13),])

dfL10 <- data.frame(parms, value = exp(l10))
dfL11 <- data.frame(parms, value = exp(l11))
dfL12 <- data.frame(parms, value = exp(l12))
dfL13 <- data.frame(parms, value = exp(l12))

#ヘシアン＝フィッシャー情報量（行列）の逆行列の符号反転
H11 <- optimHess(optp11, llmixnorm, y=y11) 
dfs1 <- sample_norm(optp11, H11, 10000)

H12 <- optimHess(optp12, llmixnorm, y=y12)
dfs2 <- sample_norm(optp12, H12, 10000)

H13 <- optimHess(optp13, llmixnorm, y=y13)
dfs3 <- sample_norm(optp13, H13, 10000)

ggplot(dfL10,aes(x=b,y=a))+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  geom_point(data=optp10, shape=4, size=1, stroke=2)+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)
ggsave("./Desktop/p_10.png")

ggplot(dfL11,aes(x=b,y=a))+
  geom_point(data=dfs1, alpha=0.1, colour="grey")+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)
ggsave("./Desktop/p_11.png")

ggplot(dfL12,aes(x=b,y=a))+
  geom_point(data=dfs2, alpha=0.1, colour="grey")+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)

ggsave("./Desktop/p_12.png")

ggplot(dfL13,aes(x=b,y=a))+
  geom_point(data=dfs3, alpha=0.1, colour="grey")+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)

ggsave("./Desktop/p_13.png")

#### Example 2
set.seed(2)
y21 <- c(rnorm(N/2,-0.4), rnorm(N/2,0.4))
x <- seq(-1, 1, length.out=N)
parms <- expand.grid(a=seq(-5,5,length.out = 200),
                     b=seq(-5,5,length.out = 200))

qplot(x, y0, geom = "line")+
  theme_bw(14)+labs(y="y (1)")
ggsave("./Desktop/p_d1.png")

ggplot()+
  geom_line(data=NULL, aes(x=x, y=y21))+
  stat_function(fun = function(x)0.4*tanh(5*x), colour="royalblue", linetype=2, size=1)+
  theme_bw(14)+labs(y="y (2)")
ggsave("./Desktop/p_d2.png")

L20 <- apply(parms, 1, lltanh, y=y0, x=x)
dfL20 <- data.frame(parms, value=exp(L20))
optp20 <- parms[which.max(L20),]

L2 <- apply(parms, 1, lltanh, y=y21, x=x)
dfL2 <- data.frame(parms, value=exp(L2))
optp2 <- unlist(parms[which.max(L2),])


H1 <- optimHess(optp2, lltanh, y=y21, x=x)
dfs <- sample_norm(optp2, H1, 10000)

ggplot(dfL20,aes(x=b,y=a))+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  geom_point(data=optp20, shape=4, size=2, stroke=1)+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)+ylim(c(-5,5))

ggsave("./Desktop/p_20.png")

ggplot(dfL2,aes(x=b,y=a))+
  geom_point(data=dfs, alpha=0.1, colour="grey")+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)
ggsave("./Desktop/p_21.png")

####Example 3
set.seed(3)
y_g <- rgamma(100,2)

parms <- expand.grid(a=seq(0.1,4,length.out = 200),
                     b=seq(0.1,4,length.out = 200))

l3 <- apply(parms, 1, llnorm, y=y_g)

optp3 <- unlist(parms[which.max(l3),])
dfL3 <- data.frame(parms, value = exp(l3))
optp3

H3 <- optimHess(optp3, llnorm, y=y_g)
dfs <- sample_norm(optp3, H3, 10000)

ggplot(dfL3,aes(x=b,y=a))+
  geom_point(data=dfs, alpha=0.1, colour="grey")+
  geom_contour(aes(z=value, colour=after_stat(level)))+
  scale_colour_binned(type = "viridis")+
  theme_bw(14)

ggsave("./Desktop/p_gamma.png")

参考文献

ベイズ統計の理論と方法

作者:渡辺澄夫
コロナ社

Amazon

In All Likelihood: Statistical Modelling And Inference Using Likelihood

作者:Pawitan, Yudi
Clarendon Press

Amazon

最尤推定の基礎

変な形の尤度関数をおもしろがるにはある程度最尤推定についての基礎知識がいるので補足する.

とはいえ, まったくの初学者が以下を理解するのことは難しいと思う. 数学的な道具としては中心極限定理とテイラー展開くらいを使う.

また, 以下での $\int f(x) \, dx$ は不定積分でなく関数の定義域全区間での積分を意味する.

さて, データ $y$ を生成したなんらかの確率分布 $q(y)$ があるとして, それを母集団と呼ぶ.

それぞれ独立に母集団 $q(y)$ に従う確率変数 $y_i$ ( $i=1,2,\ldots, n$ ) を大きさ $n$ のランダム標本（ランダムサンプル）と呼ぶ.

（現実の多くの場合には, 推測したいターゲットがこの仮想的な「母集団」と対応するように考えたりとか, ランダム標本になるようにデータのとり方を工夫したりとかするだろう.）

統計学の重要な目的なひとつは $q(y)$ におおよそ近い, 適当な確率分布 $p(y)$ をつくることである.

分析者が設定した確率分布 $p(y)$ を単にモデルと呼ぶ. パラメータ $\theta$ をつけて $p(y;\theta)$ を考え, なるべくいい $\theta$ を選ぶことにしよう.

ふたつの確率分布の近さを測るものとして, カルバック・ライブラ距離がある.

$q(y)$ と $p(y|\theta)$ のカルバック・ライブラ距離を次式で定義する.

$\displaystyle KL(p|q) = \int \log {q(y)/p(y)} q(y) \, dy\\ \displaystyle = -E_q [ \log p(y) ] + E_q[ \log q(y) ]$

ここで $E_q [ f(y)]$ は分布 $q(y)$ による期待値を表す（ $E_q [ f(y)]= \int f(y) q(y) \, dy$ ）.

カルバック・ライブラ距離の第2項が $q(y)$ のみによって定まり, モデルに依存しないので, 第1項をなるべく小さくすることを考える.

得られるのは $q(y)$ データ（確率変数） $y_i$ なので期待値を標本平均で置き換え, 対数尤度関数を次式で定義する.

$\displaystyle l_n(\theta) = \frac{1}{n} \sum_{i=1}^n \log p(y_i|\theta)$

ここではモデルに自由パラメータ $\theta$ をつけて, $p(y) =p(y|\theta)$ としている.

これはカルバック・ライブラ距離の第1項の符号反転（-1倍）なので, 対数尤度関数を最大にする $\theta$ を選ぶことで, モデルが母集団に近づくのではないかと考えている.

準備として, スコア関数 $S(\theta)$ とフィッシャー情報量 $J(\theta)$ をそれぞれ,

スコア関数：

$\displaystyle S(\theta) = \int \left(\frac{d}{d \theta} \log p(y | \theta)\right) q(y) \, dy$

フィッシャー情報量：

$\displaystyle J(\theta)=-\int \left(\frac{d^2}{d \theta^2} \log p(y | \theta)\right) q(y) \, dy$

と定義する.

次により, $S(\theta)$ の $q(y)$ による平均は 0 である.

$\displaystyle \int \left(\frac{d}{d\theta} \log p(y | \theta) \right) q(y)\, dy \\ \displaystyle = \int (\frac{d}{d\theta}p(x | \theta)) \frac{q(x)}{p(y | \theta)} \, dy \\ \displaystyle = \frac{d}{d\theta} \int q(y) \, dy\\ \displaystyle = 0.$

ここから, $S(\theta)$ の分散 $V[S(\theta)]$ は,

$\displaystyle V[S(\theta)] = E_q[{S(\theta)}^2]-E_q [{S(\theta)}]^2 = E_q[{S(\theta)}^2]$

である. 一方フィッシャー情報量は,

$\displaystyle J(\theta)= -\int \frac{d}{d \theta} \left(\left(\frac{\frac{d}{d\theta} p(y | \theta)}{p(y|\theta)}\right) q(y) \right) \, dy \\ \displaystyle =-\int \left(\frac{\frac{d^2}{d\theta^2} p(y| \theta)}{p(y|\theta)} - \frac{ ( \frac{d}{d\theta} p(x_i | \theta))^2}{p(x|\theta)^2 }\right) q(y) \, dy \\ \displaystyle =-\int \left(\frac{\frac{d^2}{d\theta^2} p(y | \theta)}{p(y|\theta)}\right) p(y|\theta) \, dx + \int \left( \frac{ ( \frac{d}{d\theta} p(y_i | \theta))^2}{p(x|\theta)^2 }\right) q(y) \, dy$

より $S(\theta)$ のときと同様, 第1項は消える. 第2項は

$\displaystyle \int \left(\frac{d}{d\theta} \log p(x | \theta)\right)^2 q(y) \, dy \\ = E_q[{S(\theta)}^2] .$

である. まとめると,

$\displaystyle V[S(\theta)] = E_q[{S(\theta)}^2]=J(\theta)$

で, $S(\theta)$ の分散は $J(\theta)$ と等しい.

さて, $l_n(\theta)$ が平均的に最大値を取る点 $\theta_0$ の周りでのテイラー展開は,

$l_n(\theta)= l_n(\theta_0) +l_n' (\theta_0) (\theta-\theta_0) + \frac{1}{2}l_n''(\theta_0)(\theta-\theta_0)^2 + \cdots$

である（テイラー展開ができることは仮定）. $h/\sqrt{n}=\theta-\theta_0$ と置くと,

$\displaystyle l_n(\theta_0 + h/\sqrt{n})= l_n(\theta_0) + \frac{l_n' (\theta_0)}{\sqrt{n}}h + \frac{1}{2n} l_n'' (\theta_0)h^2 + \cdots \\ \displaystyle \approx \frac{1}{\sqrt{n}} \left[ \sum_{i=1}^n \left( \frac{d}{d \theta} \log p(y_i | \theta) |_{\theta=\theta_0} \right)\right] h+ \frac{1}{2n}\left[ \sum_{i=1}^n \left(\frac{d^2}{d \theta^2} \log p(y_i | \theta)|_{\theta=\theta_0}\right)\right] h^2.$

上述のスコア関数とフィッシャー情報量についての性質を思い出すと, 中心極限定理より, $n$ が十分大きいとき, 標準正規分布に従う確率変数 $z$ を用いて,

$\displaystyle l_n(\theta_0 + h/\sqrt{n}) \approx z \sqrt{J(\theta_0)} h - \frac{J(\theta_0) h^2}{2} \\ \displaystyle = -\frac{J(\theta_0)}{2}\left(h - \frac{z}{\sqrt{J(\theta_0)}}\right)^2 + \frac{z^2}{2}$

という近似が成り立つ（テイラー展開の以降の項は $1/n$ より速く 0 に近づく）. $\theta=\theta_0 + h/\sqrt{n}$ であったので,

$\displaystyle l_n(\theta)= -\frac{J(\theta_0)}{2}\left(\sqrt{n}(\theta - \theta_0 )- \frac{z}{\sqrt{J(\theta_0)}}\right)^2 + \frac{z^2}{2}$

右辺は2次関数であり, $\sqrt{n}(\theta - \theta_0) = z/\sqrt{I(\theta_0)}$ のとき最大になる.

よって, $\hat \theta = \operatorname{arg\, max}_\theta l_n(\theta)$ とすると, $n$ が十分大きいとき, $\sqrt{n}(\hat \theta - \theta_0)$ は平均 0, 分散 $J(\theta_0)^{-1}$ の正規分布に従う. これが最尤法の基礎である.