混合正規分布モデルの尤度函数¶
黒木玄
2018-03-13, 2018-09-28
混合正規分布モデル
p(y|a,b)=1√2π[ae−y2/2+(1−a)e−(y−b)2/2]の尤度函数をプロットしてみよう(0≤a≤1). この確率モデルは a=1 と b=0 の場合に同一の確率分布(標準正規分布)を与える. ゆえに b=0 と a=1 の近くで通常なら起こらないことが起こることが期待される.
最尤法の漸近論によれば, 正則モデルのケースにはサンプルサイズ n→∞ で尤度函数の台はパラメーターの真の値の近くに集中して来る. しかし, 以下のプロットを見ればわかるように, パラメーターの真の値が a=0.5, b=0.1 のとき, n=218=262144 まで n を大きくしてもそのようにはならない.
この事実は, たとえ特異モデルになるパラメーターの真の値全体(今の場合は a=1 または b=0)が測度零集合であったとしても, その影響を統計学的に無視できないことを意味している. サンプルサイズ $
In [1]:
using Random: seed!
endof(a) = lastindex(a)
linspace(start, stop, length) = range(start, stop=stop, length=length)
const srand = seed!
using Distributions
using PyPlot
using Optim
d(a,b) = MixtureModel([Normal(0,1), Normal(b,1)], [a, 1-a])
#loglik(a, b, Y) = sum(logpdf(d(a,b), Y[i]) for i in 1:endof(Y)) # very slow
lpdf(a, b, y) = log(a+(1-a)*exp(b*y-b^2/2)) - y^2/2 - log(2π)/2
loglik(a, b, Y) = sum(lpdf(a, b, Y[i]) for i in 1:endof(Y))
function plot_lik(a₀, b₀, n; seed = 4649, bmin=-1.0, bmax=1.0)
srand(seed)
Y = rand(d(a₀, b₀), n)
L = 201
a = linspace(0, 1, L)
b = linspace(bmin, bmax, L)
f(a, b) = loglik(a, b, Y)
z = f.(a', b)
zmax, k = findmax(z)
#i, j = (k - 1) ÷ L + 1, mod1(k, L) # for Julia v0.6
j, i = k.I
z .= exp.(z .- zmax)
sleep(0.1)
figure(figsize=(5,5))
pcolormesh(a, b, z, cmap="CMRmap")
scatter([a₀], [b₀], label="true", color="cyan")
scatter([a[i]], [b[j]], label="MLE", color="red")
legend()
xlim(0,1)
xlabel("a")
ylabel("b")
title("\$a_0 = $a₀, b_0 = $b₀, n = $n\$")
end
Out[1]:
plot_lik (generic function with 1 method)
パラメーターの真の値が a=0.5, b=1.0 の場合¶
以下のようにサンプルサイズ n を大きくして行くと, 漸近論の予言通りに, 比較的容易に, 尤度函数の台はパラメーターの真の値の近くに集中して来る.
In [2]:
@time plot_lik(0.5, 1.0, 2^4, bmin=0, bmax=2)
2.111680 seconds (4.76 M allocations: 234.329 MiB, 4.95% gc time)
Out[2]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 1.0, n = 16$')
In [3]:
@time plot_lik(0.5, 1.0, 2^6, bmin=0, bmax=2)
0.257296 seconds (83.75 k allocations: 2.929 MiB, 2.81% gc time)
Out[3]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 1.0, n = 64$')
In [4]:
@time plot_lik(0.5, 1.0, 2^8, bmin=0, bmax=2)
0.545055 seconds (84.13 k allocations: 2.938 MiB)
Out[4]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 1.0, n = 256$')
In [5]:
@time plot_lik(0.5, 1.0, 2^10, bmin=0, bmax=2)
1.999134 seconds (86.23 k allocations: 2.993 MiB)
Out[5]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 1.0, n = 1024$')
真のパラメーターの値が a=0.5, b=0.1 の場合¶
この場合には, サンプルサイズ n を大きくしても, 漸近論の予言の通りに, 尤度函数の台がパラメーターの真の値の周囲になかなか集中して来ない.
In [6]:
@time plot_lik(0.5, 0.1, 2^4)
0.377021 seconds (661.19 k allocations: 31.127 MiB, 7.63% gc time)
Out[6]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 16$')
In [7]:
@time plot_lik(0.5, 0.1, 2^6)
0.221568 seconds (81.80 k allocations: 2.814 MiB)
Out[7]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 64$')
In [8]:
@time plot_lik(0.5, 0.1, 2^8)
0.585956 seconds (82.20 k allocations: 2.824 MiB, 12.69% gc time)
Out[8]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 256$')
In [9]:
@time plot_lik(0.5, 0.1, 2^10)
1.771081 seconds (84.23 k allocations: 2.873 MiB)
Out[9]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 1024$')
In [10]:
@time plot_lik(0.5, 0.1, 2^12)
7.434339 seconds (95.39 k allocations: 3.198 MiB)
Out[10]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 4096$')
In [11]:
@time plot_lik(0.5, 0.1, 2^14)
32.142182 seconds (132.26 k allocations: 4.042 MiB)
Out[11]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 16384$')
In [12]:
@time plot_lik(0.5, 0.1, 2^16)
118.886001 seconds (279.71 k allocations: 7.417 MiB)
Out[12]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 65536$')
In [13]:
@time plot_lik(0.5, 0.1, 2^18)
467.475714 seconds (869.61 k allocations: 20.920 MiB, 0.00% gc time)
Out[13]:
PyObject Text(0.5,1,'$a_0 = 0.5, b_0 = 0.1, n = 262144$')
In [ ]: