Instantly share code, notes, and snippets.

# genkuroki/混合正規分布について正規分布モデルで信頼区間を求めた場合.ipynb

Last active Aug 19, 2019

 { "cells": [ { "metadata": {}, "cell_type": "markdown", "source": "# 混合正規分布について正規分布モデルで信頼区間を求めた場合\n\n黒木玄\n\n2019-08-17～2019-08-19\n\nSee also https://twitter.com/genkuroki/status/1162625170034446336" }, { "metadata": { "toc": true }, "cell_type": "markdown", "source": "

\n" }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "using Distributions\nusing Plots\npyplot(size=(400,250), fmt=:svg)\nusing Random: seed!\nmyseed = 2019", "execution_count": 1, "outputs": [ { "output_type": "execute_result", "execution_count": 1, "data": { "text/plain": "2019" }, "metadata": {} } ] }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "function sim_mean(;true_dist=Normal(), α=0.05, n=100, N_iters=10^6, μ=mean(true_dist), σ=std(true_dist))\n @show true_dist\n @show μ\n @show σ\n @show α\n @show n\n @show N_iters\n M = zeros(N_iters) # sample means\n U = zeros(N_iters) # unbiased variances\n N_tdist = 0\n N_normal = 0\n z_tdist = cquantile(TDist(n-1), α/2)\n z_normal = cquantile(Normal(), α/2)\n @show z_tdist\n @show z_normal\n for i in 1:N_iters\n X = rand(true_dist, n)\n M[i] = mean(X)\n U[i] = std(X)\n N_tdist += (M[i] - z_tdist*U[i]/√n ≤ μ ≤ M[i] + z_tdist*U[i]/√n)\n N_normal += (M[i] - z_normal*σ/√n ≤ μ ≤ M[i] + z_normal*σ/√n)\n end\n P_tdist = N_tdist/N_iters\n P_normal = N_normal/N_iters\n @show P_tdist\n @show P_normal\n M, U, P_tdist, P_normal\nend", "execution_count": 2, "outputs": [ { "output_type": "execute_result", "execution_count": 2, "data": { "text/plain": "sim_mean (generic function with 1 method)" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "**sim_mean() 函数の出力の説明**\n\n* true_dist = サイズnのサンプルを生成する確率分布 = 真の分布 (デフォルト値: Normal() = 標準正規分布)\n* μ = 真の分布の平均 = 真の平均\n* σ = 真の分布の標準偏差\n* α = 有意水準 (1 - α が信頼区間の信頼係数になる) (デフォルト値: 0.05)\n* n = サンプルサイズ (デフォルト値: 100)\n* N_iters = サンプル生成を繰り返す回数 (デフォルト値: 10^6)\n* z_tdist = t分布において P(X > z) = α/2 となる z\n* z_normal = 標準正規分布において P(X > z) = α/2 となる z\n* P_tdist = t分布を使って計算した平均の信用区間に真の平均が含まれる割合\n* P_normal = 真の分布と同じ分散を持つ正規分布モデルを使って計算した平均の信用区間に真の平均が含まれる割合\n\nP_tdist と P_normal は 1 - α に近くあって欲しい.\n\nP_normal の計算では「真の分布の分散が既知」という非現実的な仮定を使っていることに注意せよ. 「真の分布の分散が未知」の場合の P_tdist との比較のために P_normal も同時に計算することにした." }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "function plotCIs(true_dist, M, U, α, n, indices=1:500)\n μ = mean(true_dist)\n z_tdist = cquantile(TDist(n-1), α/2)\n plot(; size=(1000, 150), legend=false)\n for i in indices\n A = M[i] - z_tdist*U[i]/√n\n B = M[i] + z_tdist*U[i]/√n\n color = A ≤ μ ≤ B ? :cyan : :red \n plot!([i,i], [A,B], color=color)\n end\n plot!([minimum(indices), maximum(indices)], [μ, μ], color=:black, ls=:dot)\nend\n\nfunction plot1000CIs(true_dist, M, U, α, n)\n P1 = plotCIs(true_dist, M, U, α, n, 1:500)\n P2 = plotCIs(true_dist, M, U, α, n, 501:1000)\n plot(P1, P2, size=(1000, 300), layout=grid(2,1))\nend", "execution_count": 3, "outputs": [ { "output_type": "execute_result", "execution_count": 3, "data": { "text/plain": "plot1000CIs (generic function with 1 method)" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "## 標準正規分布が真の分布の場合" }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "seed!(myseed)\nn = 100\ntrue_dist = Normal()\nM, U, P_tdist, P_normal = sim_mean(true_dist=true_dist, n=n)\nplot(x->pdf(true_dist, x), -4, 4; label=\"$true_dist\", ylim=(0, 0.5))", "execution_count": 4, "outputs": [ { "output_type": "stream", "text": "true_dist = Normal{Float64}(μ=0.0, σ=1.0)\nμ = 0.0\nσ = 1.0\nα = 0.05\nn = 100\nN_iters = 1000000\nz_tdist = 1.9842169515864176\nz_normal = 1.9599639845400592\nP_tdist = 0.949843\nP_normal = 0.949918\n", "name": "stdout" }, { "output_type": "execute_result", "execution_count": 4, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "真の分布が標準正規分布のときには, P_tdist も P_normal もほぼ理論値の 0.95 になっている. シミュレーションに伴う誤差があるのでぴったり0.95にはならない." }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "plot1000CIs(true_dist, M, U, 0.05, n)", "execution_count": 5, "outputs": [ { "output_type": "execute_result", "execution_count": 5, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "μ = mean(true_dist)\nσ = std(true_dist)\nhistogram(M; normed=true, lw=0, bins=60, label=\"sample means\", alpha=0.7)\nplot!(x->pdf(Normal(μ, σ/√n), x), label=\"normal approx\", ls=:dash)", "execution_count": 6, "outputs": [ { "output_type": "execute_result", "execution_count": 6, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "正規分布の再生性より, サンプル平均も正規分布に従う." }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "histogram((n-1)*U.^2/σ^2, normed=true, lw=0, bins=40, label=\"(n-1)U²/σ²\", alpha=0.7)\nplot!(x->pdf(Chisq(n-1), x); label=\"Chisq(n-1)\", ls=:dash)", "execution_count": 7, "outputs": [ { "output_type": "execute_result", "execution_count": 7, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "真の分布が正規分布のとき, 真の分布の分散を$\\sigma^2$と書き, サイズ$n$サンプルの不偏分散を$U^2$と書くとき,$(n-1)U^2/\\sigma^2$は自由度$n-1$のカイ二乗分布に従う." }, { "metadata": {}, "cell_type": "markdown", "source": "## 外れ値を除けば正規分布に見える混合正規分布が真の分布の場合\n\n標準正規分布に「外れ値」に見える成分を混ぜた混合正規分布の場合にはどうなるだろうか?\n\nサンプルを生成する真の分布として, 標準正規分布と平均50分散100の正規分布を$95\\% : 5\\%\$ の比で混ぜた混合正規分布を採用しよう." }, { "metadata": {}, "cell_type": "markdown", "source": "### サンプルサイズ100の場合" }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "seed!(myseed)\nn = 100\nmixnormal = MixtureModel([Normal(), Normal(50.0, 1.0)], [0.95, 0.05])\ntrue_dist = mixnormal\n\nPP = []\nfor i in 1:16\n X = rand(true_dist, n)\n P = histogram(X; label=\"\", bins=50)\n push!(PP, P)\nend\nplot(PP...; lw=0, layout=grid(4,4), size=(800, 500))", "execution_count": 8, "outputs": [ { "output_type": "execute_result", "execution_count": 8, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "@time M, U, P_tdist, P_normal = sim_mean(true_dist=true_dist, n=n)\nplot(x->pdf(true_dist, x), -10, 60; label=\"true distribution\")", "execution_count": 9, "outputs": [ { "output_type": "stream", "text": "true_dist = MixtureModel{Normal{Float64}}(K = 2)\ncomponents[1] (prior = 0.9500): Normal{Float64}(μ=0.0, σ=1.0)\ncomponents[2] (prior = 0.0500): Normal{Float64}(μ=50.0, σ=1.0)\n\nμ = 2.5\nσ = 10.943034314119645\nα = 0.05\nn = 100\nN_iters = 1000000\nz_tdist = 1.9842169515864176\nz_normal = 1.9599639845400592\nP_tdist = 0.892093\nP_normal = 0.961151\n 14.657608 seconds (214.56 M allocations: 5.982 GiB, 7.79% gc time)\n", "name": "stdout" }, { "output_type": "execute_result", "execution_count": 9, "data": { "image/svg+xml": "\r\n\r\n\r\n\r\n" }, "metadata": {} } ] }, { "metadata": {}, "cell_type": "markdown", "source": "t分布を使って機械的に計算した平均の信用区間に真の平均が含まれる割合 P_tdist が89%程度に下がってしまっていることに注意せよ. (P_normal が95%に近い値なってくれているのは, 真の分布の分散が既知であるという非現実的な仮定を使っているからである.)\n\nこのような場合があるので, サンプルの様子をプロットする場合には信用区間だけではなく, サンプルの分布そのものがある程度分かるようにしてもらえると「無駄に疑う手間」が不要になるし, サンプルの分布の様子を確認することを忘れるという初歩的な誤りも防ぐことにも役に立つだろう." }, { "metadata": { "trusted": true }, "cell_type": "code", "source": "plot1000CIs(true_dist, M, U, 0.05, n)", "execution_count": 10, "outputs": [ { "output_type": "execute_result", "execution_count": 10, "data": {
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.