Skip to content

Instantly share code, notes, and snippets.

@genkuroki

genkuroki/Animation of Bayesian estimation (Mixture Gamma Sample Case) empty.ipynb

Created November 19, 2017 08:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save genkuroki/1590ec636388f2ad9ccbbd87f53fb7e1 to your computer and use it in GitHub Desktop.
Save genkuroki/1590ec636388f2ad9ccbbd87f53fb7e1 to your computer and use it in GitHub Desktop.
Download ZIP
ベイズ推定のアニメーション (混合ガンマ分布のサンプルの場合)
Display the source blob
Display the rendered blob
Raw
Animation of Bayesian estimation (Mixture Gamma Sample Case) empty.ipynb
{
"cells": [
{
"metadata": {
"toc": "true"
},
"cell_type": "markdown",
"source": "# Table of Contents\n <p><div class=\"lev1 toc-item\"><a href=\"#ベイズ推定のアニメーション-(混合ガンマ分布のサンプルの場合)\" data-toc-modified-id=\"ベイズ推定のアニメーション-(混合ガンマ分布のサンプルの場合)-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>ベイズ推定のアニメーション (混合ガンマ分布のサンプルの場合)</a></div><div class=\"lev2 toc-item\"><a href=\"#パッケージの読み込みと諸定義\" data-toc-modified-id=\"パッケージの読み込みと諸定義-11\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>パッケージの読み込みと諸定義</a></div><div class=\"lev3 toc-item\"><a href=\"#MambaパッケージでMCMCを実行するための函数\" data-toc-modified-id=\"MambaパッケージでMCMCを実行するための函数-111\"><span class=\"toc-item-num\">1.1.1&nbsp;&nbsp;</span>MambaパッケージでMCMCを実行するための函数</a></div><div class=\"lev3 toc-item\"><a href=\"#Mambaによるシミュレーション結果のまとめの表示\" data-toc-modified-id=\"Mambaによるシミュレーション結果のまとめの表示-112\"><span class=\"toc-item-num\">1.1.2&nbsp;&nbsp;</span>Mambaによるシミュレーション結果のまとめの表示</a></div><div class=\"lev3 toc-item\"><a href=\"#WAICなどを計算するための函数\" data-toc-modified-id=\"WAICなどを計算するための函数-113\"><span class=\"toc-item-num\">1.1.3&nbsp;&nbsp;</span>WAICなどを計算するための函数</a></div><div class=\"lev3 toc-item\"><a href=\"#情報をまとめて表示するための函数\" data-toc-modified-id=\"情報をまとめて表示するための函数-114\"><span class=\"toc-item-num\">1.1.4&nbsp;&nbsp;</span>情報をまとめて表示するための函数</a></div><div class=\"lev3 toc-item\"><a href=\"#単純な正規分布モデルを最尤法で解くための函数\" data-toc-modified-id=\"単純な正規分布モデルを最尤法で解くための函数-115\"><span class=\"toc-item-num\">1.1.5&nbsp;&nbsp;</span>単純な正規分布モデルを最尤法で解くための函数</a></div><div class=\"lev2 toc-item\"><a href=\"#プロットのための函数\" data-toc-modified-id=\"プロットのための函数-12\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>プロットのための函数</a></div><div class=\"lev2 toc-item\"><a href=\"#プロットの実行-(混合ガンマ分布のサンプルの場合)\" data-toc-modified-id=\"プロットの実行-(混合ガンマ分布のサンプルの場合)-13\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>プロットの実行 (混合ガンマ分布のサンプルの場合)</a></div>"
},
{
"metadata": {},
"cell_type": "markdown",
"source": "# ベイズ推定のアニメーション (混合ガンマ分布のサンプルの場合)\n\n黒木玄\n\n2017-11-19"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "versioninfo()",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "## パッケージの読み込みと諸定義"
},
{
"metadata": {
"scrolled": false,
"trusted": true
},
"cell_type": "code",
"source": "@time begin\n #using Mamba\n \n using KernelDensity\n function makefunc_pdfkde(X)\n local ik = InterpKDE(kde(X))\n local pdfkde(x) = pdf(ik, x)\n return pdfkde\n end\n function makefunc_pdfkde(X,Y)\n local ik = InterpKDE(kde((X,Y)))\n local pdfkde(x, y) = pdf(ik, x, y)\n return pdfkde\n end\n\n using Optim\n optim_options = Optim.Options(\n store_trace = true,\n extended_trace = true\n )\n \n using QuadGK\n\n import PyPlot\n plt = PyPlot\n\n using Distributions\n @everywhere GTDist(μ, ρ, ν) = LocationScale(Float64(μ), Float64(ρ), TDist(Float64(ν)))\n @everywhere GTDist(ρ, ν) = GTDist(zero(ρ), ρ, ν)\n\n using JLD2\n using FileIO\nend\n\nusing PyCall\n@pyimport matplotlib.animation as anim\nfunction showgif(filename)\n open(filename) do f\n base64_video = base64encode(f)\n display(\"text/html\", \"\"\"<img src=\"data:image/gif;base64,$base64_video\">\"\"\")\n end\nend\n\nmacro sum(f_k, k, itr)\n quote\n begin\n local s = zero(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))\n for $(esc(k)) in $(esc(itr))\n s += $(esc(f_k))\n end\n s\n end\n end\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### MambaパッケージでMCMCを実行するための函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@everywhere mixnormal(a,b,c) = MixtureModel(Normal[Normal(b, 1.0), Normal(c, 1.0)], [1.0-a, a])\n\nmixnormal(w::AbstractVector) = mixnormal(w[1], w[2], w[3])\nunlink_mixnormal(w) = [logit(w[1]), w[2], w[3]] # unlink_mixnormal : (0,1)×R^2 -> R^3\nlink_mixnormal(z) = [invlogit(z[1]), z[2], z[3]] # link_mixnormal : R^3 → (0,1)×R^2\n\nfunction sample2model_mixnormal(Y;\n dist_model = mixnormal,\n a0 = 0.5,\n b0 = 0.0,\n c0 = 0.0,\n prior_a = Uniform(0.0, 1.0),\n prior_b = Normal(0.0, 1.0),\n prior_c = Normal(0.0, 1.0),\n chains = 2\n )\n local data = Dict{Symbol, Any}(\n :Y => Y,\n :n => length(Y),\n :a0 => a0,\n :b0 => b0,\n :c0 => c0,\n )\n \n local model = Model(\n y = Stochastic(1, (a, b, c) -> dist_model(a, b, c), false),\n a = Stochastic(() -> prior_a, true),\n b = Stochastic(() -> prior_b, true),\n c = Stochastic(() -> prior_b, true),\n )\n local scheme = [\n NUTS([:a, :b, :c])\n ]\n setsamplers!(model, scheme)\n \n local inits = [\n Dict{Symbol, Any}(\n :y => data[:Y],\n :a => data[:a0],\n :b => data[:b0],\n :c => data[:c0],\n )\n for k in 1:chains\n ]\n return model, data, inits\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@everywhere normal(mu, sigma) = Normal(mu, sigma)\n\nnormal(w::AbstractVector) = normal(w[1], w[2])\nunlink_normal(w) = [w[1], log(w[2])] # unlink_normal : R×(0,∞) -> R^2\nlink_normal(z) = [z[2], exp(z[2])] # link_normal : R^2 → R×(0,∞)\n\nfunction sample2model_normal(Y;\n dist_model = normal,\n mu0 = 0.0,\n sigma0 = 1.0,\n prior_mu = Normal(0.0, 1.0),\n prior_sigma = Truncated(Normal(1.0, 1.0), 0, Inf),\n chains = 2\n )\n local data = Dict{Symbol, Any}(\n :Y => Y,\n :n => length(Y),\n :mu0 => mu0,\n :sigma0 => sigma0,\n )\n \n local model = Model(\n y = Stochastic(1, (mu, sigma) -> dist_model(mu, sigma), false),\n mu = Stochastic(() -> prior_mu, true),\n sigma = Stochastic(() -> prior_sigma, true),\n )\n local scheme = [\n NUTS([:mu, :sigma])\n ]\n setsamplers!(model, scheme)\n \n local inits = [\n Dict{Symbol, Any}(\n :y => data[:Y],\n :mu => data[:mu0],\n :sigma => data[:sigma0],\n )\n for k in 1:chains\n ]\n return model, data, inits\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "@everywhere normal1(mu::Real) = Normal(mu, 1.0)\n\nnormal1(w::AbstractVector) = normal1(w[1])\nunlink_normal1(w) = w\nlink_normal1(z) = z\n\nfunction sample2model_normal1(Y;\n dist_model = normal1,\n mu0 = 0.0,\n prior_mu = Normal(0.0, 1.0),\n chains = 2\n )\n local data = Dict{Symbol, Any}(\n :Y => Y,\n :n => length(Y),\n :mu0 => mu0,\n )\n \n local model = Model(\n y = Stochastic(1, mu -> dist_model(mu), false),\n mu = Stochastic(() -> prior_mu, true),\n )\n local scheme = [\n NUTS([:mu])\n ]\n setsamplers!(model, scheme)\n \n local inits = [\n Dict{Symbol, Any}(\n :y => data[:Y],\n :mu => data[:mu0],\n )\n for k in 1:chains\n ]\n return model, data, inits\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### Mambaによるシミュレーション結果のまとめの表示"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "## Summary\nfunction showsummary(sim;\n sortkeys=true, figsize_t=(8, 3), figsize_c=(8, 3.5),\n show_describe=true, show_gelman=true, plot_trace=true, plot_contour=true)\n ## Summary of MCMC\n if show_describe\n println(\"\\n========== Summary:\\n\")\n display(describe(sim))\n end\n\n # Convergence Diagnostics\n if show_gelman && length(sim.chains) ≥ 2 \n println(\"========== Gelman Diagnostic:\")\n show(gelmandiag(sim))\n end\n\n ## Plot\n sleep(0.1)\n if plot_trace\n #draw(plot(sim), fmt=:png, width=10inch, height=3.5inch, nrow=1, ncol=2, ask=false)\n pyplot_trace(sim, sortkeys=sortkeys, figsize=figsize_t)\n end\n if plot_contour\n #draw(plot(sim, :contour), fmt=:png, width=10inch, height=4.5inch, nrow=1, ncol=2, ask=false)\n pyplot_contour(sim, sortkeys=sortkeys, figsize=figsize_c)\n end\nend\n\n## plot traces\nfunction pyplot_trace(sim; sortkeys = false, figsize = (8, 3))\n if sortkeys\n keys_sim = sort(keys(sim))\n else\n keys_sim = keys(sim)\n end\n for var in keys_sim\n plt.figure(figsize=figsize)\n \n plt.subplot(1,2,1)\n for k in sim.chains\n plt.plot(sim.range, sim[:,var,:].value[:,1,k], label=\"$k\", lw=0.4, alpha=0.8)\n end\n plt.xlabel(\"iterations\")\n plt.grid(ls=\":\")\n #plt.legend(loc=\"upper right\")\n plt.title(\"trace of $var\", fontsize=10)\n \n plt.subplot(1,2,2)\n local xmin = quantile(vec(sim[:,var,:].value), 0.005)\n local xmax = quantile(vec(sim[:,var,:].value), 0.995)\n for k in sim.chains\n local chain = sim[:,var,:].value[:,1,k]\n local pdfkde = makefunc_pdfkde(chain)\n local x = linspace(xmin, xmax, 201)\n plt.plot(x, pdfkde.(x), label=\"$k\", lw=0.8, alpha=0.8)\n end\n plt.xlabel(\"$var\")\n plt.grid(ls=\":\")\n plt.title(\"empirical posterior pdf of $var\", fontsize=10)\n\n plt.tight_layout()\n end\nend\n\n# plot contours\nfunction pyplot_contour(sim; sortkeys = true, figsize = (8, 3.5))\n if sortkeys\n keys_sim = sort(keys(sim))\n else\n keys_sim = keys(sim)\n end\n local m = 0\n local K = length(keys_sim)\n for i in 1:K\n for j in i+1:K\n m += 1\n mod1(m,2) == 1 && plt.figure(figsize=figsize)\n plt.subplot(1,2, mod1(m,2))\n\n local u = keys_sim[i]\n local v = keys_sim[j]\n local X = vec(sim[:,u,:].value)\n local Y = vec(sim[:,v,:].value)\n local pdfkde = makefunc_pdfkde(X,Y)\n local xmin = quantile(X, 0.005)\n local xmax = quantile(X, 0.995)\n local ymin = quantile(Y, 0.005)\n local ymax = quantile(Y, 0.995)\n local x = linspace(xmin, xmax, 200)\n local y = linspace(ymin, ymax, 200)\n \n plt.pcolormesh(x', y, pdfkde.(x',y), cmap=\"CMRmap\")\n plt.colorbar()\n plt.grid(ls=\":\")\n plt.xlabel(u)\n plt.ylabel(v)\n plt.title(\"posterior of ($u, $v)\", fontsize=10)\n\n mod1(m,2) == 2 && plt.tight_layout()\n\n if 2*m == K*(K-1) && mod1(m,2) == 1\n plt.subplot(1,2,2)\n \n plt.pcolormesh(y', x, pdfkde.(x,y'), cmap=\"CMRmap\")\n plt.colorbar()\n plt.grid(ls=\":\")\n plt.xlabel(v)\n plt.ylabel(u)\n plt.title(\"posterior of ($v, $u)\", fontsize=10)\n\n plt.tight_layout()\n end\n end\n end\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### WAICなどを計算するための函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# loglik[l,i] = lpdf(w_l, Y_i) と chain[l,:] = w_l を取り出す函数を作る函数\n#\nfunction makefunc_loglikchainof(lpdf, symbols, Y)\n #\n # loglikchainof(sim) で loglik[l,i] と chain[l,:] が抽出される\n #\n local function loglikchainof(sim)\n local val = sim[:, symbols, :].value\n local chain = vcat((val[:,:,k] for k in 1:size(val,3))...)\n local L = size(chain,1)\n local n = length(Y)\n local loglik = Array{Float64, 2}(L, n)\n for i in 1:n\n for l in 1:L\n loglik[l,i] = lpdf(chain[l,:], Y[i])\n end\n end\n return loglik, chain\n end\n return loglikchainof\nend\n\n# 予測分布函数 p^*(x,y) = mean of { lpdf(w_l, y) }_{l=1}^L を作る函数\n#\nfunction makefunc_pdfpred(lpdf, chain)\n local L = size(chain,1)\n local pred_Bayes(y) = @sum(exp(lpdf((@view chain[l,:]), y)), l, 1:L)/L\n return pred_Bayes\nend\n\n# loglik[l,i] からWAICを計算する函数\n#\nfunction WAICof(loglik)\n local L, n\n L, n = size(loglik)\n local T_n = -mean(log(mean(exp(loglik[l,i]) for l in 1:L)) for i in 1:n)\n local V_n = sum(var(loglik[:,i], corrected=false) for i in 1:n)\n local WAIC = 2*n*T_n + 2*V_n\n return WAIC, 2*n*T_n, 2*V_n\nend\n\n# loglik[l,i] からLOOCVを素朴に計算する函数\n#\nfunction LOOCVof(loglik)\n local L, n\n L, n = size(loglik)\n local LOOCV = 2*sum(log(mean(exp(-loglik[l,i]) for l in 1:L)) for i in 1:n)\n return LOOCV\nend\n\n# 自由エネルギー(の2倍)を計算するための函数\n# \n# 自由エネルギーの2の逆温度微分は E^β_w[2n L_n] に等しいので、\n# それを β=0.0 から 1.0 まで数値積分すれば自由エネルギーの2倍を計算できる。\n#\nfunction FreeEnergyof(loglik)\n local E2nLn = makefunc_E2nLn(loglik)\n local F = quadgk(E2nLn, 0.0, 1.0)[1]\n return F, E2nLn\nend\n\nfunction makefunc_E2nLn(loglik)\n local L = size(loglik)[1]\n local negloglik = -sum(loglik, 2)\n local negloglik_n = negloglik .- maximum(negloglik)\n local function E2nLn(beta)\n local Edenominator = @sum( exp((1-beta)*negloglik_n[l]), l, 1:L)/L\n if Edenominator == zero(Edenominator) || !isfinite(Edenominator)\n return zero(Edenominator)\n end\n local Enumerator = @sum(negloglik[l]*exp((1-beta)*negloglik_n[l]), l, 1:L)/L\n return 2*Enumerator/Edenominator\n end\n return E2nLn\nend\n\n# loglik[l,i] からWBICを計算する函数\n#\nfunction WBICof(loglik)\n local E2nLn = makefunc_E2nLn(loglik)\n local n = size(loglik, 2)\n local WBIC = E2nLn(1/log(n))\n return WBIC\nend\n\n# 汎化損失を計算する函数\n#\nfunction GeneralizationLossof(pdftrue, pdfpred; xmin=-10.0, xmax=10.0)\n local f(x) = -pdftrue(x)*log(pdfpred(x))\n return quadgk(f, xmin, xmax)\nend\n\n# Shannon情報量を計算する函数\n#\nShannonInformationof(pdftrue; xmin=-10.0, xmax=10.0) = GeneralizationLossof(pdftrue, pdftrue, xmin=xmin, xmax=xmax)\nShannonInformationof(dist::Distribution, n) = 2*n*entropy(dist)\n\n# 最尤法を実行して AIC と BIC を計算する函数\n#\n# lpdf(w,y) = log p(y|w)\n#\n# w = link_model(z) ←実ベクトル z 全体をパラメーター w の空間内に移す函数\n# z = unlink_model(w) ←link_model(z)の逆函数\n# これらは optimize() 函数の適用時にパラメーター w が定義域から外れることを防ぐための函数達\n#\n# (X[i], Y[i] はサンプル\n#\n# chain は loglik, chain = loglikchainof(sim) で作った chain\n#\n# optimize函数は1変数函数の場合には初期条件を与えて解を求める函数ではなくなるので、\n# その場合には2変数函数に拡張してから使用している.\n#\nfunction AICandBICof(lpdf, link_model, unlink_model, Y, chain)\n local n = length(Y)\n local L = size(chain,1)\n local nparams = size(chain,2)\n local negloglik(z) = -@sum(lpdf(link_model(z), Y[i]), i, 1:n)\n local minnegloglik_chain, l\n minnegloglik_chain, l = findmin(negloglik(unlink_model(chain[l,:])) for l in 1:L)\n local o, minnegloglik, param_AIC\n if size(chain,2) == 1\n local f(z) = negloglik([z[1]]) + z[2]^2/eps()\n o = optimize(f, [unlink_model(chain[l,:])[1], 0.0])\n minnegloglik = o.minimum\n param_AIC = link_model([o.minimizer[1]])\n else\n o = optimize(negloglik, unlink_model(chain[l,:]))\n minnegloglik = o.minimum\n param_AIC = link_model(o.minimizer)\n end\n local T_AIC = 2.0*minnegloglik\n local V_AIC = 2.0*nparams\n local T_BIC = T_AIC\n local V_BIC = nparams*log(n)\n local AIC = T_AIC + V_AIC\n local BIC = T_BIC + V_BIC\n return AIC, T_AIC, V_AIC, \n BIC, T_BIC, V_BIC, \n param_AIC\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 情報をまとめて表示するための函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function statsof(sim, Y; \n dist_true=mixnormal(0.0, 0.0, 0.0), \n dist_model=mixnormal, link_model=link_mixnormal, unlink_model=unlink_mixnormal)\n \n local n = length(Y)\n \n local lpdf(w, y) = logpdf(dist_model(w), y)\n \n local loglikchainof = makefunc_loglikchainof(lpdf, sort(keys(sim)), Y)\n local loglik, chain \n loglik, chain = loglikchainof(sim)\n \n local WAIC, T_WAIC, V_WAIC\n WAIC, T_WAIC, V_WAIC = WAICof(loglik)\n local LOOCV = LOOCVof(loglik)\n \n local WBIC = WBICof(loglik)\n local FreeEnergy = FreeEnergyof(loglik)[1]\n \n local param_Bayes = vec(mean(chain, 1))\n local pred_Bayes = makefunc_pdfpred(lpdf, chain)\n \n local GeneralizationLoss = 2*n*GeneralizationLossof(x->pdf(dist_true,x), pred_Bayes)[1]\n \n local AIC, T_AIC, V_AIC, BIC, T_BIC, V_BIC, param_MLE\n AIC, T_AIC, V_AIC, BIC, T_BIC, V_BIC, param_MLE = AICandBICof(lpdf, link_model, unlink_model, Y, chain)\n \n local pred_MLE(y) = exp(lpdf(param_MLE, y))\n return WAIC, T_WAIC, V_WAIC, LOOCV, GeneralizationLoss,\n WBIC, FreeEnergy,\n param_Bayes, pred_Bayes, \n AIC, T_AIC, V_AIC, \n BIC, T_BIC, V_BIC,\n param_MLE, pred_MLE\nend\n\nfunction show_all_results(dist_true, Y, sim; statsfunc=statsof, \n dist_model=mixnormal, link_model=link_mixnormal, unlink_model=link_mixnormal,\n figsize=(6,4.2), xmin=-4.0, xmax=6.0)\n WAIC, T_WAIC, V_WAIC, LOOCV, GeneralizationLoss,\n WBIC, FreeEnergy,\n param_Bayes, pred_Bayes, \n AIC, T_AIC, V_AIC, \n BIC, T_BIC, V_BIC,\n param_MLE, pred_MLE = statsfunc(sim, Y, dist_true=dist_true, \n dist_model=dist_model, link_model=link_model, unlink_model=unlink_model)\n \n n = length(Y)\n println(\"\\n=== Estimates by $dist_model (n = $n) ===\")\n @show param_Bayes\n @show param_MLE\n \n println(\"--- Information Criterions\")\n println(\"* AIC = $AIC = $T_AIC + $V_AIC\")\n println(\"* GenLoss = $GeneralizationLoss\")\n println(\"* WAIC = $WAIC = $T_WAIC + $V_WAIC\")\n println(\"* LOOCV = $LOOCV\")\n println(\"---\")\n println(\"* BIC = $BIC = $T_BIC + $V_BIC\")\n println(\"* FreeEnergy = $FreeEnergy\")\n println(\"* WBIC = $WBIC\")\n\n println(\"=\"^78 * \"\\n\")\n \n sleep(0.1)\n plt.figure(figsize=figsize)\n kde_sample = makefunc_pdfkde(Y)\n x = linspace(xmin, xmax, 201)\n plt.plot(x, pdf.(dist_true, x), label=\"true distribution\")\n plt.scatter(Y, kde_sample.(Y), label=\"sample\", s=10, color=\"k\", alpha=0.5)\n plt.plot(x, kde_sample.(x), label=\"KDE of sample\", color=\"k\", alpha=0.5)\n plt.plot(x, pred_Bayes.(x), label=\"Baysian predictive\", ls=\"--\")\n plt.plot(x, pred_MLE.(x), label=\"MLE predictive\", ls=\"-.\")\n plt.xlabel(\"x\")\n plt.ylabel(\"probability density\")\n plt.grid(ls=\":\")\n plt.legend(fontsize=8)\n plt.title(\"Estimates by $dist_model: n = $(length(Y))\")\nend\n\nfunction plotsample(dist_true, Y; figsize=(6,4.2), xmin=-4.0, xmax=6.0)\n sleep(0.1)\n plt.figure(figsize=figsize)\n kde_sample = makefunc_pdfkde(Y)\n x = linspace(xmin, xmax, 201)\n plt.plot(x, pdf.(dist_true, x), label=\"true dist.\")\n plt.scatter(Y, kde_sample.(Y), label=\"sample\", s=10, color=\"k\", alpha=0.5)\n plt.plot(x, kde_sample.(x), label=\"KDE of sample\", color=\"k\", alpha=0.5)\n plt.xlabel(\"x\")\n plt.ylabel(\"probability density\")\n plt.grid(ls=\":\")\n plt.legend(fontsize=8)\n plt.title(\"Sample size n = $(length(Y))\")\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "### 単純な正規分布モデルを最尤法で解くための函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# dist_true 分布に従う乱数で生成したサンプルの配列 Y を与えると\n# 正規分布モデルの最尤法で推定して、AICなどを返す函数\n#\nfunction fit_Normal(dist_true, Y)\n local n = length(Y)\n local d = fit(Normal,Y)\n local mu = d.μ\n local sigma = d.σ\n local pred_Normal(y) = pdf(d,y)\n local T_Normal = -2*sum(logpdf.(d,Y))\n local V_Normal = 4.0\n local AIC_Normal = T_Normal + V_Normal\n local TB_Normal = T_Normal\n local VB_Normal = 2.0*log(n)\n local BIC_Normal = TB_Normal + VB_Normal\n local f(y) = -pdf(dist_true, y)*logpdf(d, y)\n local GL_Normal = 2*n*quadgk(f, -10, 10)[1]\n return mu, sigma, pred_Normal, \n AIC_Normal, T_Normal, V_Normal, \n BIC_Normal, TB_Normal, VB_Normal, \n GL_Normal\nend\n\n# グラフをプロットする範囲が xmin から xmax まで\n#\nfunction show_fit_Normal(dist_true, Y; figsize=(6,4.2), xmin=-4.0, xmax=6.0)\n mu, sigma, pred_Normal, \n AIC_Normal, T_Normal, V_Normal, \n BIC_Normal, TB_Normal, VB_Normal, \n GL_Normal = fit_Normal(dist_true, Y)\n println(\"--- Normal Fitting\")\n println(\"* μ = $mu\")\n println(\"* σ = $sigma\")\n println(\"* GenLoss = $GL_Normal\")\n println(\"* AIC = $AIC_Normal = $T_Normal + $V_Normal\")\n println(\"* BIC = $BIC_Normal = $TB_Normal + $VB_Normal\")\n \n sleep(0.1)\n plt.figure(figsize=figsize)\n kde_sample = makefunc_pdfkde(Y)\n x = linspace(xmin, xmax, 201)\n plt.plot(x, pdf.(dist_true, x), label=\"true distribution\")\n plt.scatter(Y, kde_sample.(Y), label=\"sample\", s=10, color=\"k\", alpha=0.5)\n plt.plot(x, kde_sample.(x), label=\"KDE of sample\", color=\"k\", alpha=0.5)\n plt.plot(x, pred_Normal.(x), label=\"Normal predictive\", ls=\"--\")\n plt.xlabel(\"x\")\n plt.ylabel(\"probability density\")\n plt.grid(ls=\":\")\n plt.legend(fontsize=8)\n plt.title(\"Sample size n = $(length(Y))\")\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function InformationCriterions(dist_true, Y, sim; dist_model=mixnormal)\n local n = length(Y)\n local lpdf(w, y) = logpdf(dist_model(w), y)\n\n local loglikchainof = makefunc_loglikchainof(lpdf, sort(keys(sim)), Y)\n local loglik, chain \n loglik, chain = loglikchainof(sim)\n\n local WAIC, T_WAIC, V_WAIC\n WAIC, T_WAIC, V_WAIC = WAICof(loglik)\n local LOOCV = LOOCVof(loglik)\n \n local pred_Bayes = makefunc_pdfpred(lpdf, chain)\n local GL = 2*n*GeneralizationLossof(x->pdf(dist_true,x), pred_Bayes)[1]\n \n local WBIC = WBICof(loglik)\n local FreeEnergy = FreeEnergyof(loglik)[1]\n\n return chain, WAIC, T_WAIC, V_WAIC, LOOCV, GL, WBIC, FreeEnergy\n end",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true
},
"cell_type": "markdown",
"source": "## プロットのための函数"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "function plotsim(t; modelname=\"normal\", xmin=-3, xmax=8, y1max=0.45, y2max=0.40)\n plt.clf()\n \n local n = t\n local Y = Sample[1:n]\n local Chain = eval(Symbol(:Chain_, modelname))\n local chain = Chain[n]\n local dist_model = eval(Symbol(modelname))\n local lpdf(w,x) = logpdf(dist_model(w),x)\n local pred = makefunc_pdfpred(lpdf, chain)\n local x = linspace(xmin,xmax,201)\n local kl = eval(Symbol(:kl_, modelname))\n local wt = eval(Symbol(:wt_, modelname))\n \n plt.subplot(121)\n plt.plot(x, pdf.(dist_true,x), color=\"black\", ls=\":\", label=\"true dist\")\n plt.scatter(Y, pdf.(dist_true,Y), color=\"red\", s=10, alpha=0.5, label=\"sample\") \n plt.plot(x, pred.(x), label=\"predictive\")\n plt.ylim(-y1max/40, y1max)\n plt.grid(ls=\":\")\n plt.legend()\n plt.title(\"$modelname: n = $n\")\n \n plt.subplot(122)\n plt.plot(1:n, kl[1:n], label=\"KL information\")\n plt.plot(1:n, wt[1:n], label=\"(WAIC \\$-\\$ T\\$_\\\\mathrm{true}\\$)/2n\")\n plt.xlabel(\"sample size n\")\n plt.ylabel(\"\\$-\\$ log(probability)\")\n plt.xlim(-1, N+1)\n plt.ylim(min(minimum(kl), minimum(wt[5:end]))-y2max/40, y2max)\n plt.grid(ls=\":\")\n plt.legend()\n plt.title(\"$modelname: n = $n\")\n \n plt.tight_layout()\n plt.plot()\nend",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true
},
"cell_type": "markdown",
"source": "## プロットの実行 (混合ガンマ分布のサンプルの場合)\n\nGIFアニメーションを作成したい場合には適当にコメントアウトを外して実行する。"
},
{
"metadata": {
"collapsed": true
},
"cell_type": "markdown",
"source": "なぜか normal モデルの場合の計算がうまく行っていない.\n\nサンプルサイズが 50 をちょっと超えたあたりから, 事後予測分布がおかしな形で動かなくなる.\n\nこれを書いている時点(2017-11-19)の時点で原因不明."
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "N0 = 53",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "seed = 4649\n#seed = 12345\n#seed = 2017\n\n#dataname = \"Gamma128Sample$seed\"\ndataname = \"MixGamma128Sample$seed\"\n\ndata = load(\"$dataname.jld2\")\nfor v in sort(collect(keys(data)))\n ex = parse(\"$v = data[\\\"$v\\\"]\")\n # println(ex)\n eval(ex)\nend\n\n@show mu_true = mean(dist_true)\n@show sigma_true = std(dist_true)\n@show dist_normal_fitting = Normal(mu_true, sigma_true)\n\nshannon_true = Shannon[1]/2\ngl_normal_fitting = GeneralizationLossof(x->pdf(dist_true,x), x->pdf(dist_normal_fitting,x))[1]\nkl_normal_fitting = gl_normal_fitting - shannon_true\n\nKL_normal = GL_normal - Shannon\nKL_normal1 = GL_normal1 - Shannon\nKL_mixnormal = GL_mixnormal - Shannon\n\nWT_normal = WAIC_normal - T_true\nWT_normal1 = WAIC_normal1 - T_true\nWT_mixnormal = WAIC_mixnormal - T_true\n\nkl_normal = [KL_normal[n]/(2n) for n in 1:N]\nkl_normal1 = [KL_normal1[n]/(2n) for n in 1:N]\nkl_mixnormal = [KL_mixnormal[n]/(2n) for n in 1:N]\n\nwt_normal = [WT_normal[n]/(2n) for n in 1:N]\nwt_normal1 = [WT_normal1[n]/(2n) for n in 1:N]\nwt_mixnormal = [WT_mixnormal[n]/(2n) for n in 1:N]\n\n@show shannon_true\n@show gl_normal_fitting\n@show kl_normal_fitting\n@show exp(-kl_normal_fitting)\nprintln()\n@show kl_normal\n@show kl_normal1\n@show kl_mixnormal\n\nsleep(0.1)\n\nN, Ntmp = N0, N\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal\", xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nN = Ntmp\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal1\", xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"mixnormal\", xmin=-1, xmax=11, y1max=0.35)\n\nN, Ntmp = N0, N\nmodelname = \"normal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\nN = Ntmp\n\nmodelname = \"normal1\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\n\nmodelname = \"mixnormal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "seed = 4649\n#seed = 12345\n#seed = 2017\n\n#dataname = \"Gamma128Sample$seed\"\ndataname = \"MixGamma128Sample$seed\"\n\ndata = load(\"$dataname.jld2\")\nfor v in sort(collect(keys(data)))\n ex = parse(\"$v = data[\\\"$v\\\"]\")\n # println(ex)\n eval(ex)\nend\n\n@show mu_true = mean(dist_true)\n@show sigma_true = std(dist_true)\n@show dist_normal_fitting = Normal(mu_true, sigma_true)\n\nshannon_true = Shannon[1]/2\ngl_normal_fitting = GeneralizationLossof(x->pdf(dist_true,x), x->pdf(dist_normal_fitting,x))[1]\nkl_normal_fitting = gl_normal_fitting - shannon_true\n\nKL_normal = GL_normal - Shannon\nKL_normal1 = GL_normal1 - Shannon\nKL_mixnormal = GL_mixnormal - Shannon\n\nWT_normal = WAIC_normal - T_true\nWT_normal1 = WAIC_normal1 - T_true\nWT_mixnormal = WAIC_mixnormal - T_true\n\nkl_normal = [KL_normal[n]/(2n) for n in 1:N]\nkl_normal1 = [KL_normal1[n]/(2n) for n in 1:N]\nkl_mixnormal = [KL_mixnormal[n]/(2n) for n in 1:N]\n\nwt_normal = [WT_normal[n]/(2n) for n in 1:N]\nwt_normal1 = [WT_normal1[n]/(2n) for n in 1:N]\nwt_mixnormal = [WT_mixnormal[n]/(2n) for n in 1:N]\n\n@show shannon_true\n@show gl_normal_fitting\n@show kl_normal_fitting\n@show exp(-kl_normal_fitting)\nprintln()\n@show kl_normal\n@show kl_normal1\n@show kl_mixnormal\n\nsleep(0.1)\n\nN, Ntmp = N0, N\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal\", xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nN = Ntmp\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal1\", xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"mixnormal\", xmin=-1, xmax=11, y1max=0.35)\n\nN, Ntmp = N0, N\nmodelname = \"normal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\nN = Ntmp\n\nmodelname = \"normal1\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\n\nmodelname = \"mixnormal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "#seed = 4649\nseed = 12345\n#seed = 2017\n\n#dataname = \"Gamma128Sample$seed\"\ndataname = \"MixGamma128Sample$seed\"\n\ndata = load(\"$dataname.jld2\")\nfor v in sort(collect(keys(data)))\n ex = parse(\"$v = data[\\\"$v\\\"]\")\n # println(ex)\n eval(ex)\nend\n\n@show mu_true = mean(dist_true)\n@show sigma_true = std(dist_true)\n@show dist_normal_fitting = Normal(mu_true, sigma_true)\n\nshannon_true = Shannon[1]/2\ngl_normal_fitting = GeneralizationLossof(x->pdf(dist_true,x), x->pdf(dist_normal_fitting,x))[1]\nkl_normal_fitting = gl_normal_fitting - shannon_true\n\nKL_normal = GL_normal - Shannon\nKL_normal1 = GL_normal1 - Shannon\nKL_mixnormal = GL_mixnormal - Shannon\n\nWT_normal = WAIC_normal - T_true\nWT_normal1 = WAIC_normal1 - T_true\nWT_mixnormal = WAIC_mixnormal - T_true\n\nkl_normal = [KL_normal[n]/(2n) for n in 1:N]\nkl_normal1 = [KL_normal1[n]/(2n) for n in 1:N]\nkl_mixnormal = [KL_mixnormal[n]/(2n) for n in 1:N]\n\nwt_normal = [WT_normal[n]/(2n) for n in 1:N]\nwt_normal1 = [WT_normal1[n]/(2n) for n in 1:N]\nwt_mixnormal = [WT_mixnormal[n]/(2n) for n in 1:N]\n\n@show shannon_true\n@show gl_normal_fitting\n@show kl_normal_fitting\n@show exp(-kl_normal_fitting)\nprintln()\n@show kl_normal\n@show kl_normal1\n@show kl_mixnormal\n\nsleep(0.1)\n\nN, Ntmp = N0, N\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal\", xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nN = Ntmp\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal1\", xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"mixnormal\", xmin=-1, xmax=11, y1max=0.35)\n\nN, Ntmp = N0, N\nmodelname = \"normal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\nN = Ntmp\n\nmodelname = \"normal1\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\n\nmodelname = \"mixnormal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "#seed = 4649\n#seed = 12345\nseed = 2017\n\n#dataname = \"Gamma128Sample$seed\"\ndataname = \"MixGamma128Sample$seed\"\n\ndata = load(\"$dataname.jld2\")\nfor v in sort(collect(keys(data)))\n ex = parse(\"$v = data[\\\"$v\\\"]\")\n # println(ex)\n eval(ex)\nend\n\n@show mu_true = mean(dist_true)\n@show sigma_true = std(dist_true)\n@show dist_normal_fitting = Normal(mu_true, sigma_true)\n\nshannon_true = Shannon[1]/2\ngl_normal_fitting = GeneralizationLossof(x->pdf(dist_true,x), x->pdf(dist_normal_fitting,x))[1]\nkl_normal_fitting = gl_normal_fitting - shannon_true\n\nKL_normal = GL_normal - Shannon\nKL_normal1 = GL_normal1 - Shannon\nKL_mixnormal = GL_mixnormal - Shannon\n\nWT_normal = WAIC_normal - T_true\nWT_normal1 = WAIC_normal1 - T_true\nWT_mixnormal = WAIC_mixnormal - T_true\n\nkl_normal = [KL_normal[n]/(2n) for n in 1:N]\nkl_normal1 = [KL_normal1[n]/(2n) for n in 1:N]\nkl_mixnormal = [KL_mixnormal[n]/(2n) for n in 1:N]\n\nwt_normal = [WT_normal[n]/(2n) for n in 1:N]\nwt_normal1 = [WT_normal1[n]/(2n) for n in 1:N]\nwt_mixnormal = [WT_mixnormal[n]/(2n) for n in 1:N]\n\n@show shannon_true\n@show gl_normal_fitting\n@show kl_normal_fitting\n@show exp(-kl_normal_fitting)\nprintln()\n@show kl_normal\n@show kl_normal1\n@show kl_mixnormal\n\nsleep(0.1)\n\nN, Ntmp = N0, N\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal\", xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nN = Ntmp\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"normal1\", xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\n\nplt.figure(figsize=(10,3.5))\nplotsim(N, modelname=\"mixnormal\", xmin=-1, xmax=11, y1max=0.35)\n\nN, Ntmp = N0, N\nmodelname = \"normal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\nN = Ntmp\n\nmodelname = \"normal1\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.45, y2max=1.0)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()\n\nmodelname = \"mixnormal\"\nfile = \"$dataname$modelname.gif\"\nplot1frame(t) = plotsim(t, modelname=modelname, xmin=-1, xmax=11, y1max=0.35)\nfig = plt.figure(figsize=(10,3.5))\nframes = [N; 1;1;1;1; 2;2;2;3;3;3;4;4;4; 5;5;5;6;6;7;7;8;8;9;9;10;10; 11:N; N;N;N;N;N;N;N;N;N]\ninterval = 100\n@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)\n#@time myanim[:save](file, writer=\"imagemagick\")\nsleep(0.1)\nshowgif(file)\nplt.clf()",
"execution_count": null,
"outputs": []
},
{
"metadata": {
"collapsed": true,
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"_draft": {
"nbviewer_url": "https://gist.github.com/d688c6e5e8c427dba9856d37a12bbee1"
},
"gist": {
"id": "d688c6e5e8c427dba9856d37a12bbee1",
"data": {
"description": "ベイズ推定のアニメーション (混合ガンマ分布のサンプルの場合)",
"public": true
}
},
"kernelspec": {
"name": "julia-0.6.0",
"display_name": "Julia 0.6.0",
"language": "julia"
},
"language_info": {
"file_extension": ".jl",
"name": "julia",
"mimetype": "application/julia",
"version": "0.6.0"
},
"toc": {
"threshold": 4,
"number_sections": true,
"toc_cell": true,
"toc_window_display": false,
"toc_section_display": "block",
"sideBar": true,
"navigate_menu": true,
"moveMenuLeft": true,
"widenNotebook": false,
"colors": {
"hover_highlight": "#DAA520",
"selected_highlight": "#FFD700",
"running_highlight": "#FF0000",
"wrapper_background": "#FFFFFF",
"sidebar_border": "#EEEEEE",
"navigate_text": "#333333",
"navigate_num": "#000000"
},
"nav_menu": {
"height": "12px",
"width": "252px"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment