混合正規分布モデルの尤度函数

黒木玄

2018-03-13, 2018-09-28

混合正規分布モデル

$$ p(y|a,b) = \frac{1}{\sqrt{2\pi}}\left[ ae^{-y^2/2} + (1-a)e^{-(y-b)^2/2} \right] $$

の尤度函数をプロットしてみよう($0\le a\le 1$). この確率モデルは $a=1$ と $b=0$ の場合に同一の確率分布(標準正規分布)を与える. ゆえに $b=0$ と $a=1$ の近くで通常なら起こらないことが起こることが期待される.

最尤法の漸近論によれば, 正則モデルのケースにはサンプルサイズ $n\to\infty$ で尤度函数の台はパラメーターの真の値の近くに集中して来る. しかし, 以下のプロットを見ればわかるように, パラメーターの真の値が $a=0.5$, $b=0.1$ のとき, $n = 2^{18} = 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 [ ]: