using Base.Meta: parse using BenchmarkTools using PyPlot using SymPy using SpecialFunctions using Statistics const latex = sympy[:latex] using LaTeXStrings latexstring_displaystyle(args...; kwargs...) = LaTeXString(raw"$$" * prod(latex.(args; kwargs...)) * raw"$$") latexdisp(args...; kwargs...) = display(latexstring_displaystyle(args...; kwargs...)) const ls = latexstring_displaystyle const ld = latexdisp @syms a b c d @syms x y μ μ₀ real=true Y = [symbols("Y$i", real=true) for i in 1:5] @syms σ λ λ₀ θ α β ξ η positive=true Y p_N = 1/√(2PI)*exp(-λ*(y-μ)^2/2)*√λ p_μ = p_N(μ=>μ₀, y=>μ, λ=>λ*λ₀) p_λ = exp(-λ/θ)*λ^(α-1)/(gamma(α)*θ^α) # for simplicity, put μ₀ = 0. p_μ = p_μ(μ₀=>0) logp_N = simplify(log(p_N)) logp_μ = simplify(log(p_μ)) logp_λ = expand(log(p_λ)) p_N = simplify(p_N) p_μ = simplify(p_μ) p_λ = simplify(p_λ) logp_N = expand(logp_N) logp_μ = expand(logp_μ) logp_λ = expand(logp_λ) [p_N, p_μ, p_λ, logp_N, logp_μ, logp_λ] @time I_μ = simplify(integrate(expand(logp_N*p_μ), (μ, -oo, oo))) # 対数尤度函数の期待値 @time Elogp_N = simplify(integrate(expand(I_μ*p_λ), (λ, 0, oo))) A1 = -2Elogp_N A2 = log(2PI) + α*θ*y^2 + 1/λ₀ - (polygamma(Sym(0),α) + log(θ)) simplify(A1 - A2) expand(logp_N^2) @time J_μ = integrate(expand(logp_N^2*p_μ), (μ, -oo, oo)) # 対数尤度函数の二乗の期待値 @time Elogp_N2 = simplify(integrate(expand(J_μ*p_λ), (λ, 0, oo))) @syms ψ₀ ψ₁ B1 = 4Elogp_N2(polygamma(Sym(0),α)=>ψ₀, polygamma(Sym(1),α)=>ψ₁) B2 = ( (log(2PI))^2 + (α+1)*α*θ^2*y^4 + 6α*θ/λ₀*y^2 + 3/λ₀^2 + ψ₁ + (ψ₀ + log(θ))^2 + 2log(2PI)*(α*θ*y^2 + 1/λ₀) - 2log(2PI)*(ψ₀ + log(θ)) - 2( y^2*(θ + α*θ*(ψ₀ + log(θ))) + (ψ₀+log(θ))/λ₀) ) simplify(B1 - B2) p_μ = p_μ(μ=>μ-μ₀) n = length(Y) Ybar = mean(Y) Yvar = sum((Y.-Ybar).^2)/length(Y) # 共役事前分布 prior = simplify(p_μ*p_λ) p5_N = simplify(prod(p_N(y=>v)^β for v in Y)) # 次の p5 を分配函数で割ったものが事後分布になる. p5 = simplify(p5_N*p_μ*p_λ) # 以下の式を μ, λ に関する定数倍を除いて上の共役事前分布と等しくするためには, # 共役事前分布の4つのパラメーターをどのように変えなければいけないかを以下で計算する. # n=5 の場合の対数尤度函数 logp5 = expand(log(p5)) # exp() の中の式における -λ の係数 c_λ = -coeff(collect(logp5, λ), λ) @syms a b c s = solve(c_λ - a/2*(μ-b)^2 - c, [a,b,c])[1] ld("a = ", s[a]) ld("b = ", s[b]) ld("c = ", s[c]) # μ に関する平方完成 simplify(c_λ - s[a]/2*(μ-s[b])^2 - s[c]) # α₅ は対数尤度函数における log(λ) の係数に 1/2 を足したもの α₅ = coeff(collect(logp5, log(λ)), log(λ)) + Sym(1)/2 λ₀₅ = s[a] μ₀₅ = s[b] θ₅ = 1/simplify(s[c]) α₅ - (α + n*β/2) λ₀₅ - (λ₀+n*β) simplify(μ₀₅ - (λ₀*μ₀+n*β*Ybar)/(λ₀+n*β)) tmp = (1/θ)*(1 + (n*β*θ/2)*(λ₀/(λ₀+n*β)*(Ybar - μ₀)^2 + Yvar)) simplify(1/θ₅ - tmp) prior_param = (μ₀=>0.0, θ=>1/α, λ₀=>0.001, α=>0.000001) posterior_param = (μ₀=>μ₀₅, λ₀=>λ₀₅, α=>α₅, θ=>θ₅) data_sample = (Y[1]=>0.7, Y[2]=>1.1, Y[3]=>1.2, Y[4]=>1.2, Y[5]=>1.5) prior = (p_μ*p_λ)(prior_param...) posterior = (p_μ*p_λ)(posterior_param...)(β=>1.0, prior_param..., data_sample...) posterior = posterior(β=>1.0, prior_param..., data_sample...) posterior = simplify(posterior) # lambdify版 f_prior = lambdify(prior, [μ, λ]) f_posterior = lambdify(posterior, [μ, λ]) μs = -1:0.03:2 λs = 10 .^(-1:0.04:3) z_prior = f_prior.(μs', λs) z_posterior = f_posterior.(μs', λs) sleep(0.1) figure(figsize=(10,4)) subplot(121) pcolormesh(μs, λs, z_prior, cmap="CMRmap") colorbar() grid(ls=":") xlabel("μ") ylabel("λ") yscale("log") title("prior") subplot(122) pcolormesh(μs, λs, z_posterior, cmap="CMRmap") colorbar() grid(ls=":") xlabel("μ") ylabel("λ") yscale("log") title("posterior") tight_layout() # eval parse 版 eval(parse("g_prior(μ,λ) = $prior")) eval(parse("g_posterior(μ,λ) = $posterior")) μs = -1:0.03:2 λs = 10 .^(-1:0.04:3) z_prior = g_prior.(μs', λs) z_posterior = g_posterior.(μs', λs) sleep(0.1) figure(figsize=(10,4)) subplot(121) pcolormesh(μs, λs, z_prior, cmap="CMRmap") colorbar() grid(ls=":") xlabel("μ") ylabel("λ") yscale("log") title("prior") subplot(122) pcolormesh(μs, λs, z_posterior, cmap="CMRmap") colorbar() grid(ls=":") xlabel("μ") ylabel("λ") yscale("log") title("posterior") tight_layout() # SymPy expression posterior # String version of SymPy expression "$posterior" # 101×101 ≈ 10^4 (μ, λ)'s μs = collect(-1:0.03:2) λs = 10 .^(-1:0.04:3) @show length.((μs,λs)); # SymPy lambdify 版 f_posterior = lambdify(posterior, [μ, λ]) @benchmark z_posterior = f_posterior.(μs', λs) # Julia eval parse 版 eval(parse("g_posterior(μ,λ) = $posterior")) @benchmark z_posterior = g_posterior.(μs', λs) N(posterior) # Julia eval parse N 版 eval(parse("h_posterior(μ,λ) = $(N(posterior))")) @benchmark z_posterior = h_posterior.(μs', λs) # SymPy lambdify 版 meshgrid ケース μ_grid = repeat(μs, 1, length(λs))' λ_grid = repeat(λs, 1, length(μs)) f_posterior = lambdify(posterior, [μ, λ]) @benchmark z_posterior = f_posterior.(μ_grid, λ_grid) figure(figsize=(5,4)) pcolormesh(μs, λs, z_posterior, cmap="CMRmap") colorbar() grid(ls=":") xlabel("μ") ylabel("λ") yscale("log") title("posterior")