# 正規分布の共役事前分布関連の公式のSymPyによる導出¶

2018-04-02, 2018-12-20

In [1]:
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

Out[1]:
latexdisp (generic function with 1 method)
In [2]:
@syms a b c d
@syms x y μ μ₀ real=true
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 lambdify版と Julia eval parse版の比較¶ In [46]: # SymPy expression posterior  Out[46]: \begin{equation*}\frac{0.00953627593040395 \sqrt{2} λ^{2.000001} e^{- λ \left(2.5005 \left(μ - 0.56994300569943\right)^{2} + 0.166650670065987\right)}}{\sqrt{\pi}}\end{equation*} In [47]: # String version of SymPy expression "$posterior"

Out[47]:
"0.00953627593040395*sqrt(2)*λ^2.000001*exp(-λ*(2.5005*(μ - 0.56994300569943)^2 + 0.166650670065987))/sqrt(pi)"
In [48]:
# 101×101 ≈ 10^4 (μ, λ)'s
μs = collect(-1:0.03:2)
λs = 10 .^(-1:0.04:3)
@show length.((μs,λs));

length.((μs, λs)) = (101, 101)

In [49]:
# SymPy lambdify 版
f_posterior = lambdify(posterior, [μ, λ])
@benchmark z_posterior = f_posterior.(μs', λs)

Out[49]:
BenchmarkTools.Trial:
memory estimate:  1.01 MiB
allocs estimate:  51017
--------------
minimum time:     2.796 ms (0.00% GC)
median time:      2.912 ms (0.00% GC)
mean time:        3.189 ms (4.28% GC)
maximum time:     64.073 ms (92.81% GC)
--------------
samples:          1564
evals/sample:     1
In [50]:
# Julia eval parse 版
eval(parse("g_posterior(μ,λ) = $posterior")) @benchmark z_posterior = g_posterior.(μs', λs)  Out[50]: BenchmarkTools.Trial: memory estimate: 79.89 KiB allocs estimate: 5 -------------- minimum time: 1.162 ms (0.00% GC) median time: 1.183 ms (0.00% GC) mean time: 1.241 ms (2.11% GC) maximum time: 62.258 ms (98.05% GC) -------------- samples: 4028 evals/sample: 1 In [51]: N(posterior)  Out[51]: \begin{equation*}\frac{0.00953627593040395 \sqrt{2} λ^{2.000001} e^{- λ \left(2.5005 \left(μ - 0.56994300569943\right)^{2} + 0.166650670065987\right)}}{\sqrt{\pi}}\end{equation*} In [52]: # Julia eval parse N 版 eval(parse("h_posterior(μ,λ) =$(N(posterior))"))
@benchmark z_posterior = h_posterior.(μs', λs)

Out[52]:
BenchmarkTools.Trial:
memory estimate:  79.89 KiB
allocs estimate:  5
--------------
minimum time:     1.148 ms (0.00% GC)
median time:      1.220 ms (0.00% GC)
mean time:        1.377 ms (2.50% GC)
maximum time:     82.969 ms (98.31% GC)
--------------
samples:          3621
evals/sample:     1

ちなみにJuliaにおいて meshgrid は次のセルのようにして実現できる.

しかし, meshgridは, 同一のデータを大量に複製することによって, メモリを無駄に使用する. 可能ならば避けた方がコンピューターの効率的な利用という観点からは合理的である.

In [53]:
# 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)

Out[53]:
BenchmarkTools.Trial:
memory estimate:  1.01 MiB
allocs estimate:  51016
--------------
minimum time:     2.809 ms (0.00% GC)
median time:      2.862 ms (0.00% GC)
mean time:        3.066 ms (4.38% GC)
maximum time:     65.052 ms (95.40% GC)
--------------
samples:          1630
evals/sample:     1
In [54]:
figure(figsize=(5,4))
pcolormesh(μs, λs, z_posterior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("posterior")

Out[54]:
PyObject Text(0.5,1,'posterior')

## 続く¶

かもしれない。

In [ ]: