Julia言語のSymPy.jlで変分ベイズの例題を理解する

黒木玄

2018-04-02, 2018-12-20

2018-04-01 PythonのSymPyで変分ベイズの例題を理解する (StatModeling Memorandum) と同じことを, Julia言語からPython SymPyを利用してやってみよう.

このブログ記事は私が今までに見たSymPyの使い方に関する記事の中で最高のものであった. このノートでは詳しい解説は省略するが, ブログ記事の方には各実行結果に詳細なコメントが付けられており, SymPyを使いこなしたい人にとっては必見のブログ記事だと思った.

ほとんどの節は上のブログ記事のPython版のコードのJulia言語への翻訳である. 「変数変換をして再計算」の節と最後の「変分近似無しの計算」の節は私自身による.

JuliaでSymPyを使いたい人は tutorial2016-01-22 Julia+SymPyの勝手なまとめ (たけし備忘録) を参照するとよい.

In [1]:
using Base.Meta: parse
using SpecialFunctions
linspace(a,b,L) = range(a, stop=b, length=L)

using PyPlot
using PyCall
using SymPy

# https://docs.sympy.org/latest/modules/printing.html#sympy.printing.julia.julia_code
const julia_code = sympy[:julia_code]

@show versioninfo()
println()
@show PyCall.pyprogramname
@show PyCall.pyversion
@show PyCall.conda
@show PyCall.libpython
@show sympy[:__version__];
Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, haswell)
Environment:
  JULIA_CMDSTAN_HOME = C:\CmdStan
  JULIA_NUM_THREADS = 4
  JULIA_PKGDIR = C:\JuliaPkg
versioninfo() = nothing

PyCall.pyprogramname = "C:\\Anaconda3\\python.exe"
PyCall.pyversion = v"3.6.5"
PyCall.conda = false
PyCall.libpython = "C:\\Anaconda3\\python36.dll"
sympy[:__version__] = "1.3"

SymPyの変数の準備

In [2]:
y = symbols("y", real=true)
t = symbols("t", positive=true)
μ, μ₀ = symbols("μ μ₀", real=true)
Y = symbols("Y1:4", real=true)
T = symbols("T1:4", positive=true)
τ, λ₀, a₀, b₀ = symbols("τ λ₀ a₀ b₀", positive=true)
σ, σ² = symbols("σ σ²", positive=true)
θ, α = symbols("θ α", positive=true);

確率密度函数を直接定義

In [3]:
pdf_Normal(mu, s2, t) = exp(-(t-mu)^2/(2*s2))/√(2*PI*s2)
@show pdf_Normal(μ, σ², y)
pdf_Normal(μ, σ², y) = sqrt(2)*exp(-(y - μ)^2/(2*σ²))/(2*sqrt(pi)*sqrt(σ²))
Out[3]:
\begin{equation*}\frac{\sqrt{2} e^{- \frac{\left(y - μ\right)^{2}}{2 σ²}}}{2 \sqrt{\pi} \sqrt{σ²}}\end{equation*}
In [4]:
pdf_Gamma(alpha, theta, t) = exp(-t/theta)*t^(alpha-1)/(gamma(alpha)*theta^alpha)
@show pdf_Gamma(α, θ, y)
pdf_Gamma(α, θ, y) = y^(α - 1)*θ^(-α)*exp(-y/θ)/gamma(α)
Out[4]:
\begin{equation*}\frac{y^{α - 1} θ^{- α} e^{- \frac{y}{θ}}}{\Gamma\left(α\right)}\end{equation*}
In [5]:
p_y = pdf_Normal(μ, 1/τ, y)
@show p_y
p_y = sqrt(2)*sqrt(τ)*exp(-τ*(y - μ)^2/2)/(2*sqrt(pi))
Out[5]:
\begin{equation*}\frac{\sqrt{2} \sqrt{τ} e^{- \frac{τ \left(y - μ\right)^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [6]:
p_μ = pdf_Normal(μ₀, 1/(λ₀*τ), μ)
@show p_μ
p_μ = sqrt(2)*sqrt(λ₀)*sqrt(τ)*exp(-λ₀*τ*(μ - μ₀)^2/2)/(2*sqrt(pi))
Out[6]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ₀} \sqrt{τ} e^{- \frac{λ₀ τ \left(μ - μ₀\right)^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [7]:
p_τ = pdf_Gamma(a₀, 1/b₀, τ)
@show p_τ
p_τ = b₀^a₀*τ^(a₀ - 1)*exp(-b₀*τ)/gamma(a₀)
Out[7]:
\begin{equation*}\frac{b₀^{a₀} τ^{a₀ - 1} e^{- b₀ τ}}{\Gamma\left(a₀\right)}\end{equation*}

積分計算の例

In [8]:
I1 = integrate(p_μ, (μ, -oo,  oo))
@show I1
I1 = 1
Out[8]:
\begin{equation*}1\end{equation*}
In [9]:
I2 = integrate(p_τ, (τ, 0, oo))
@show I2
I2 = b₀^a₀*b₀^(-a₀ + 1)/b₀
Out[9]:
\begin{equation*}\frac{b₀^{a₀} b₀^{- a₀ + 1}}{b₀}\end{equation*}
In [10]:
I2 = simplify(I2)
Out[10]:
\begin{equation*}1\end{equation*}

正規分布モデルの共役事前分布に関する事前予測分布

まず, $\mu$ について積分する.

In [11]:
Z1_μ = simplify(integrate(p_y*p_μ, (μ, -oo, oo)))
Out[11]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ₀} \sqrt{τ} e^{\frac{τ \left(- y^{2} + 2 y μ₀ - μ₀^{2} + \frac{\left(y - μ₀\right)^{2}}{λ₀ + 1}\right)}{2}}}{2 \sqrt{\pi} \sqrt{λ₀ + 1}}\end{equation*}
In [12]:
F1_μ = simplify(-log(Z1_μ))
Out[12]:
\begin{equation*}\frac{τ \left(- \left(y - μ₀\right)^{2} + \left(λ₀ + 1\right) \left(y^{2} - 2 y μ₀ + μ₀^{2}\right)\right)}{2 \left(λ₀ + 1\right)} - \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (τ \right )}}{2} + \frac{\log{\left (λ₀ + 1 \right )}}{2} + \frac{\log{\left (2 \right )}}{2} + \frac{\log{\left (\pi \right )}}{2}\end{equation*}

以下では $\mu_0=0$, $t=y$ と置いてから, $\tau$ に関する積分を実行する. $t$ を $y-\mu_0$ に置き換えれば $\mu_0$ を復活させることができる. SymPyにとって積分計算を易しくなるように $t>0$ と仮定してある.

In [13]:
F1_0_μ = F1_μ(μ₀=>0, y=>t)
Out[13]:
\begin{equation*}\frac{τ \left(t^{2} \left(λ₀ + 1\right) - t^{2}\right)}{2 \left(λ₀ + 1\right)} - \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (τ \right )}}{2} + \frac{\log{\left (λ₀ + 1 \right )}}{2} + \frac{\log{\left (2 \right )}}{2} + \frac{\log{\left (\pi \right )}}{2}\end{equation*}
In [14]:
F1_0_μ = expand(F1_0_μ)
Out[14]:
\begin{equation*}\frac{t^{2} λ₀ τ}{2 \left(λ₀ + 1\right)} - \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (τ \right )}}{2} + \frac{\log{\left (λ₀ + 1 \right )}}{2} + \frac{\log{\left (2 \right )}}{2} + \frac{\log{\left (\pi \right )}}{2}\end{equation*}
In [15]:
Z1_0_μ = simplify(exp(-F1_0_μ))
Out[15]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ₀} \sqrt{τ} e^{- \frac{t^{2} λ₀ τ}{2 \left(λ₀ + 1\right)}}}{2 \sqrt{\pi} \sqrt{λ₀ + 1}}\end{equation*}
In [16]:
Z1_0 = integrate(Z1_0_μ*p_τ, (τ, 0, oo))
Out[16]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ₀} \left(1 + \frac{t^{2} λ₀}{2 b₀ \left(λ₀ + 1\right)}\right)^{- a₀ - \frac{1}{2}} \Gamma\left(a₀ + \frac{1}{2}\right)}{2 \sqrt{\pi} \sqrt{b₀} \sqrt{λ₀ + 1} \Gamma\left(a₀\right)}\end{equation*}

これが正規分布モデルを共役事前分布で平均して得た事前予測分布の式である($t=y-\mu_0$). t分布になっている.

同時分布の対数 ($\log p$) の準備

In [17]:
log_p = sum(log(p_y(y=>x)) for x in Y) + log(p_μ) + log(p_τ)
log_p = simplify(log_p)
@show log_p
log_p = -Y1^2*τ/2 + Y1*μ*τ - Y2^2*τ/2 + Y2*μ*τ - Y3^2*τ/2 + Y3*μ*τ + a₀*log(b₀) + a₀*log(τ) - b₀*τ - λ₀*μ^2*τ/2 + λ₀*μ*μ₀*τ - λ₀*μ₀^2*τ/2 - 3*μ^2*τ/2 + log(a₀) + log(λ₀)/2 + log(τ) - log(gamma(a₀ + 1)) - 2*log(pi) - 2*log(2)
Out[17]:
\begin{equation*}- \frac{Y_{1}^{2} τ}{2} + Y_{1} μ τ - \frac{Y_{2}^{2} τ}{2} + Y_{2} μ τ - \frac{Y_{3}^{2} τ}{2} + Y_{3} μ τ + a₀ \log{\left (b₀ \right )} + a₀ \log{\left (τ \right )} - b₀ τ - \frac{λ₀ μ^{2} τ}{2} + λ₀ μ μ₀ τ - \frac{λ₀ μ₀^{2} τ}{2} - \frac{3 μ^{2} τ}{2} + \log{\left (a₀ \right )} + \frac{\log{\left (λ₀ \right )}}{2} + \log{\left (τ \right )} - \log{\left (\Gamma\left(a₀ + 1\right) \right )} - 2 \log{\left (\pi \right )} - 2 \log{\left (2 \right )}\end{equation*}

上の式の書き直し:

$$ \begin{aligned} &- \frac{Y_{1}^{2} τ}{2} + Y_{1} μ τ - \frac{Y_{2}^{2} τ}{2} + Y_{2} μ τ - \frac{Y_{3}^{2} τ}{2} + Y_{3} μ τ + a₀ \log{\left (b₀ \right )} + a₀ \log{\left (τ \right )} - b₀ τ \\ &- \frac{λ₀ τ}{2} μ^{2} + λ₀ μ μ₀ τ - \frac{λ₀ τ}{2} μ₀^{2} - \frac{3 τ}{2} μ^{2} + \log{\left (a₀ \right )} + \frac{1}{2} \log{\left (λ₀ \right )} + \log{\left (τ \right )} - \log{\left (\Gamma{\left(a₀ + 1 \right)} \right )} \\ &- 2 \log{\left (\pi \right )} - 2 \log{\left (2 \right )} \end{aligned} $$
In [18]:
log_p_for_μ = integrate(diff(log_p, μ), μ)
@show log_p_for_μ
log_p_for_μ = μ^2*(-λ₀*τ/2 - 3*τ/2) + μ*(Y1*τ + Y2*τ + Y3*τ + λ₀*μ₀*τ)
Out[18]:
\begin{equation*}μ^{2} \left(- \frac{λ₀ τ}{2} - \frac{3 τ}{2}\right) + μ \left(Y_{1} τ + Y_{2} τ + Y_{3} τ + λ₀ μ₀ τ\right)\end{equation*}
In [19]:
expand(log_p_for_μ)
Out[19]:
\begin{equation*}Y_{1} μ τ + Y_{2} μ τ + Y_{3} μ τ - \frac{λ₀ μ^{2} τ}{2} + λ₀ μ μ₀ τ - \frac{3 μ^{2} τ}{2}\end{equation*}
In [20]:
collect(expand(log_p_for_μ), μ)
Out[20]:
\begin{equation*}μ^{2} \left(- \frac{λ₀ τ}{2} - \frac{3 τ}{2}\right) + μ \left(Y_{1} τ + Y_{2} τ + Y_{3} τ + λ₀ μ₀ τ\right)\end{equation*}
In [21]:
log_p_for_τ = integrate(diff(log_p, τ), τ)
@show log_p_for_τ
log_p_for_τ = τ*(-Y1^2 + 2*Y1*μ - Y2^2 + 2*Y2*μ - Y3^2 + 2*Y3*μ - 2*b₀ - λ₀*μ^2 + 2*λ₀*μ*μ₀ - λ₀*μ₀^2 - 3*μ^2)/2 + (a₀ + 1)*log(τ)
Out[21]:
\begin{equation*}\frac{τ \left(- Y_{1}^{2} + 2 Y_{1} μ - Y_{2}^{2} + 2 Y_{2} μ - Y_{3}^{2} + 2 Y_{3} μ - 2 b₀ - λ₀ μ^{2} + 2 λ₀ μ μ₀ - λ₀ μ₀^{2} - 3 μ^{2}\right)}{2} + \left(a₀ + 1\right) \log{\left (τ \right )}\end{equation*}
In [22]:
expand(log_p_for_τ)
Out[22]:
\begin{equation*}- \frac{Y_{1}^{2} τ}{2} + Y_{1} μ τ - \frac{Y_{2}^{2} τ}{2} + Y_{2} μ τ - \frac{Y_{3}^{2} τ}{2} + Y_{3} μ τ + a₀ \log{\left (τ \right )} - b₀ τ - \frac{λ₀ μ^{2} τ}{2} + λ₀ μ μ₀ τ - \frac{λ₀ μ₀^{2} τ}{2} - \frac{3 μ^{2} τ}{2} + \log{\left (τ \right )}\end{equation*}
In [23]:
log_p_for_τ = collect(expand(log_p_for_τ), τ)
Out[23]:
\begin{equation*}a₀ \log{\left (τ \right )} + τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ - \frac{Y_{2}^{2}}{2} + Y_{2} μ - \frac{Y_{3}^{2}}{2} + Y_{3} μ - b₀ - \frac{λ₀ μ^{2}}{2} + λ₀ μ μ₀ - \frac{λ₀ μ₀^{2}}{2} - \frac{3 μ^{2}}{2}\right) + \log{\left (τ \right )}\end{equation*}

できるところまで解析的に求める

q1_μ

In [24]:
log_q1_μ = integrate(log_p_for_μ * p_τ, (τ, 0, oo))
@show log_q1_μ
log_q1_μ = Y1*μ*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y2*μ*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y3*μ*gamma(a₀ + 1)/(b₀*gamma(a₀)) - λ₀*μ^2*gamma(a₀ + 1)/(2*b₀*gamma(a₀)) + λ₀*μ*μ₀*gamma(a₀ + 1)/(b₀*gamma(a₀)) - 3*μ^2*gamma(a₀ + 1)/(2*b₀*gamma(a₀))
Out[24]:
\begin{equation*}\frac{Y_{1} μ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{2} μ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{3} μ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} - \frac{λ₀ μ^{2} \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)} + \frac{λ₀ μ μ₀ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} - \frac{3 μ^{2} \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)}\end{equation*}
In [25]:
log_q1_μ = collect(expand(log_q1_μ), μ)
@show log_q1_μ
log_q1_μ = μ^2*(-λ₀*gamma(a₀ + 1)/(2*b₀*gamma(a₀)) - 3*gamma(a₀ + 1)/(2*b₀*gamma(a₀))) + μ*(Y1*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y2*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y3*gamma(a₀ + 1)/(b₀*gamma(a₀)) + λ₀*μ₀*gamma(a₀ + 1)/(b₀*gamma(a₀)))
Out[25]:
\begin{equation*}μ^{2} \left(- \frac{λ₀ \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)} - \frac{3 \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)}\right) + μ \left(\frac{Y_{1} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{2} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{3} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{λ₀ μ₀ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)}\right)\end{equation*}
In [26]:
q1_μ = exp(log_q1_μ)
z1_μ = simplify(integrate(q1_μ, (μ, -oo, oo)))
@show z1_μ
z1_μ = sqrt(2)*sqrt(pi)*sqrt(b₀)*exp(Y1^2*a₀/(2*b₀*(λ₀ + 3)) + Y1*Y2*a₀/(b₀*(λ₀ + 3)) + Y1*Y3*a₀/(b₀*(λ₀ + 3)) + Y1*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) + Y2^2*a₀/(2*b₀*(λ₀ + 3)) + Y2*Y3*a₀/(b₀*(λ₀ + 3)) + Y2*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) + Y3^2*a₀/(2*b₀*(λ₀ + 3)) + Y3*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) + a₀*λ₀^2*μ₀^2/(2*b₀*(λ₀ + 3)))/(sqrt(a₀)*sqrt(λ₀ + 3))
Out[26]:
\begin{equation*}\frac{\sqrt{2} \sqrt{\pi} \sqrt{b₀} e^{\frac{Y_{1}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} + \frac{Y_{1} Y_{2} a₀}{b₀ \left(λ₀ + 3\right)} + \frac{Y_{1} Y_{3} a₀}{b₀ \left(λ₀ + 3\right)} + \frac{Y_{1} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} + \frac{Y_{2}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} + \frac{Y_{2} Y_{3} a₀}{b₀ \left(λ₀ + 3\right)} + \frac{Y_{2} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} + \frac{Y_{3}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} + \frac{Y_{3} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} + \frac{a₀ λ₀^{2} μ₀^{2}}{2 b₀ \left(λ₀ + 3\right)}}}{\sqrt{a₀} \sqrt{λ₀ + 3}}\end{equation*}
In [27]:
q1_μ = 1/z1_μ * exp(log_q1_μ)
@show q1_μ
q1_μ = sqrt(2)*sqrt(a₀)*sqrt(λ₀ + 3)*exp(μ^2*(-λ₀*gamma(a₀ + 1)/(2*b₀*gamma(a₀)) - 3*gamma(a₀ + 1)/(2*b₀*gamma(a₀))) + μ*(Y1*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y2*gamma(a₀ + 1)/(b₀*gamma(a₀)) + Y3*gamma(a₀ + 1)/(b₀*gamma(a₀)) + λ₀*μ₀*gamma(a₀ + 1)/(b₀*gamma(a₀))))*exp(-Y1^2*a₀/(2*b₀*(λ₀ + 3)) - Y1*Y2*a₀/(b₀*(λ₀ + 3)) - Y1*Y3*a₀/(b₀*(λ₀ + 3)) - Y1*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) - Y2^2*a₀/(2*b₀*(λ₀ + 3)) - Y2*Y3*a₀/(b₀*(λ₀ + 3)) - Y2*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) - Y3^2*a₀/(2*b₀*(λ₀ + 3)) - Y3*a₀*λ₀*μ₀/(b₀*(λ₀ + 3)) - a₀*λ₀^2*μ₀^2/(2*b₀*(λ₀ + 3)))/(2*sqrt(pi)*sqrt(b₀))
Out[27]:
\begin{equation*}\frac{\sqrt{2} \sqrt{a₀} \sqrt{λ₀ + 3} e^{μ^{2} \left(- \frac{λ₀ \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)} - \frac{3 \Gamma\left(a₀ + 1\right)}{2 b₀ \Gamma\left(a₀\right)}\right) + μ \left(\frac{Y_{1} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{2} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{Y_{3} \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)} + \frac{λ₀ μ₀ \Gamma\left(a₀ + 1\right)}{b₀ \Gamma\left(a₀\right)}\right)} e^{- \frac{Y_{1}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} - \frac{Y_{1} Y_{2} a₀}{b₀ \left(λ₀ + 3\right)} - \frac{Y_{1} Y_{3} a₀}{b₀ \left(λ₀ + 3\right)} - \frac{Y_{1} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} - \frac{Y_{2}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} - \frac{Y_{2} Y_{3} a₀}{b₀ \left(λ₀ + 3\right)} - \frac{Y_{2} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} - \frac{Y_{3}^{2} a₀}{2 b₀ \left(λ₀ + 3\right)} - \frac{Y_{3} a₀ λ₀ μ₀}{b₀ \left(λ₀ + 3\right)} - \frac{a₀ λ₀^{2} μ₀^{2}}{2 b₀ \left(λ₀ + 3\right)}}}{2 \sqrt{\pi} \sqrt{b₀}}\end{equation*}

q1_τ (何も工夫せずに計算できてしまった!しかし非常に汚い)

In [28]:
log_q1_τ = integrate(log_p_for_τ * p_μ, (μ, -oo, oo))
@show log_q1_τ
log_q1_τ = -Y1^2*τ/2 + Y1*μ₀*τ - Y2^2*τ/2 + Y2*μ₀*τ - Y3^2*τ/2 + Y3*μ₀*τ + a₀*log(τ) - b₀*τ - 3*μ₀^2*τ/2 + log(τ) - 1/2 - 3/(2*λ₀)
Out[28]:
\begin{equation*}- \frac{Y_{1}^{2} τ}{2} + Y_{1} μ₀ τ - \frac{Y_{2}^{2} τ}{2} + Y_{2} μ₀ τ - \frac{Y_{3}^{2} τ}{2} + Y_{3} μ₀ τ + a₀ \log{\left (τ \right )} - b₀ τ - \frac{3 μ₀^{2} τ}{2} + \log{\left (τ \right )} - \frac{1}{2} - \frac{3}{2 λ₀}\end{equation*}
In [29]:
log_q1_τ = integrate(diff(log_q1_τ, τ), τ)
@show log_q1_τ
log_q1_τ = τ*(-Y1^2 + 2*Y1*μ₀ - Y2^2 + 2*Y2*μ₀ - Y3^2 + 2*Y3*μ₀ - 2*b₀ - 3*μ₀^2)/2 + (a₀ + 1)*log(τ)
Out[29]:
\begin{equation*}\frac{τ \left(- Y_{1}^{2} + 2 Y_{1} μ₀ - Y_{2}^{2} + 2 Y_{2} μ₀ - Y_{3}^{2} + 2 Y_{3} μ₀ - 2 b₀ - 3 μ₀^{2}\right)}{2} + \left(a₀ + 1\right) \log{\left (τ \right )}\end{equation*}
In [30]:
log_q1_τ = collect(log_q1_τ, τ)
@show log_q1_τ
log_q1_τ = τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2) + (a₀ + 1)*log(τ)
Out[30]:
\begin{equation*}τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right) + \left(a₀ + 1\right) \log{\left (τ \right )}\end{equation*}
In [31]:
q1_τ = logcombine(exp(log_q1_τ))
@show q1_τ
q1_τ = τ^(a₀ + 1)*exp(τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2))
Out[31]:
\begin{equation*}τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\end{equation*}
In [32]:
z1_τ = integrate(q1_τ, (τ, 0, oo))
@show z1_τ
z1_τ = Piecewise((2*(exp_polar(I*pi)*polar_lift(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2))^(-a₀ - 1)*gamma(a₀ + 2)/(Y1^2 - 2*Y1*μ₀ + Y2^2 - 2*Y2*μ₀ + Y3^2 - 2*Y3*μ₀ + 2*b₀ + 3*μ₀^2), pi/2 > Abs(arg(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2) + pi)), (Integral(τ^(a₀ + 1)*exp(τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2)), (τ, 0, oo)), True))
Out[32]:
\begin{equation*}\begin{cases} \frac{2 \left(e^{i \pi} \operatorname{polar\_lift}{\left (- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2} \right )}\right)^{- a₀ - 1} \Gamma\left(a₀ + 2\right)}{Y_{1}^{2} - 2 Y_{1} μ₀ + Y_{2}^{2} - 2 Y_{2} μ₀ + Y_{3}^{2} - 2 Y_{3} μ₀ + 2 b₀ + 3 μ₀^{2}} & \text{for}\: \frac{\pi}{2} > \left|{\arg{\left (- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2} \right )} + \pi}\right| \\\int\limits_{0}^{\infty} τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\, dτ & \text{otherwise} \end{cases}\end{equation*}

↑何も工夫せずに積分がそのまま計算できてしまった! しかし非常に汚い!

In [33]:
q1_τ = 1/z1_τ * q1_τ
@show q1_τ
q1_τ = τ^(a₀ + 1)*Piecewise(((exp_polar(I*pi)*polar_lift(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2))^(a₀ + 1)*(Y1^2 - 2*Y1*μ₀ + Y2^2 - 2*Y2*μ₀ + Y3^2 - 2*Y3*μ₀ + 2*b₀ + 3*μ₀^2)/(2*gamma(a₀ + 2)), pi/2 > Abs(arg(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2) + pi)), (1/Integral(τ^(a₀ + 1)*exp(τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2)), (τ, 0, oo)), True))*exp(τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2))
Out[33]:
\begin{equation*}τ^{a₀ + 1} \left(\begin{cases} \frac{\left(e^{i \pi} \operatorname{polar\_lift}{\left (- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2} \right )}\right)^{a₀ + 1} \left(Y_{1}^{2} - 2 Y_{1} μ₀ + Y_{2}^{2} - 2 Y_{2} μ₀ + Y_{3}^{2} - 2 Y_{3} μ₀ + 2 b₀ + 3 μ₀^{2}\right)}{2 \Gamma\left(a₀ + 2\right)} & \text{for}\: \frac{\pi}{2} > \left|{\arg{\left (- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2} \right )} + \pi}\right| \\\frac{1}{\int\limits_{0}^{\infty} τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\, dτ} & \text{otherwise} \end{cases}\right) e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\end{equation*}

変数変換をして再計算

計算できたが, 出力汚な過ぎるので, $\mu_0=0$, $Y_i=T_i>0$ と置いてやり直してみる.

In [34]:
q2_τ = logcombine(exp(log_q1_τ))
@show q2_τ
q2_τ = τ^(a₀ + 1)*exp(τ*(-Y1^2/2 + Y1*μ₀ - Y2^2/2 + Y2*μ₀ - Y3^2/2 + Y3*μ₀ - b₀ - 3*μ₀^2/2))
Out[34]:
\begin{equation*}τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\end{equation*}
In [35]:
q2_0_τ = subs(q2_τ, (μ₀, 0), zip(Y,T)...)
Out[35]:
\begin{equation*}τ^{a₀ + 1} e^{τ \left(- \frac{T_{1}^{2}}{2} - \frac{T_{2}^{2}}{2} - \frac{T_{3}^{2}}{2} - b₀\right)}\end{equation*}
In [36]:
z2_0_τ = integrate(q2_0_τ, (τ, 0, oo))
@show z2_0_τ
z2_0_τ = (T3^2/(2*(T1^2/2 + T2^2/2 + b₀)) + 1)^(-a₀ - 2)*(T1^2/2 + T2^2/2 + b₀)^(-a₀ - 1)*gamma(a₀ + 2)/(T1^2/2 + T2^2/2 + b₀)
Out[36]:
\begin{equation*}\frac{\left(\frac{T_{3}^{2}}{2 \left(\frac{T_{1}^{2}}{2} + \frac{T_{2}^{2}}{2} + b₀\right)} + 1\right)^{- a₀ - 2} \left(\frac{T_{1}^{2}}{2} + \frac{T_{2}^{2}}{2} + b₀\right)^{- a₀ - 1} \Gamma\left(a₀ + 2\right)}{\frac{T_{1}^{2}}{2} + \frac{T_{2}^{2}}{2} + b₀}\end{equation*}
In [37]:
z2_0_τ = simplify(z2_0_τ)
@show z2_0_τ
z2_0_τ = 2^(a₀ + 2)*(T1^2 + T2^2 + T3^2 + 2*b₀)^(-a₀ - 2)*gamma(a₀ + 2)
Out[37]:
\begin{equation*}2^{a₀ + 2} \left(T_{1}^{2} + T_{2}^{2} + T_{3}^{2} + 2 b₀\right)^{- a₀ - 2} \Gamma\left(a₀ + 2\right)\end{equation*}
In [38]:
q2_0_τ = 1/z2_0_τ * q2_0_τ
@show q2_0_τ
q2_0_τ = 2^(-a₀ - 2)*τ^(a₀ + 1)*(T1^2 + T2^2 + T3^2 + 2*b₀)^(a₀ + 2)*exp(τ*(-T1^2/2 - T2^2/2 - T3^2/2 - b₀))/gamma(a₀ + 2)
Out[38]:
\begin{equation*}\frac{2^{- a₀ - 2} τ^{a₀ + 1} \left(T_{1}^{2} + T_{2}^{2} + T_{3}^{2} + 2 b₀\right)^{a₀ + 2} e^{τ \left(- \frac{T_{1}^{2}}{2} - \frac{T_{2}^{2}}{2} - \frac{T_{3}^{2}}{2} - b₀\right)}}{\Gamma\left(a₀ + 2\right)}\end{equation*}

q1_τ の計算に再挑戦

In [39]:
q1_τ = logcombine(exp(log_q1_τ))
display(q1_τ)
\begin{equation*}τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}\end{equation*}
In [40]:
coef = coeff(collect(log_q1_τ, τ), τ)
display(coef)

sol = solve(diff(coef, Y[1]), Y[1])[1]
display(sol) #-> μ₀

replacements = [(var, sol) for var in Y]
display(subs(coef, replacements...)) #-> -b₀
\begin{equation*}- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\end{equation*}
\begin{equation*}μ₀\end{equation*}
\begin{equation*}- b₀\end{equation*}
In [41]:
ξ = symbols("ξ", positive=true)
z1_τ = simplify(integrate(τ^(a₀+1)*exp(-ξ*τ), (τ, 0, oo)))
z1_τ = subs(z1_τ, (ξ, -coef))
@show z1_τ
z1_τ = (Y1^2/2 - Y1*μ₀ + Y2^2/2 - Y2*μ₀ + Y3^2/2 - Y3*μ₀ + b₀ + 3*μ₀^2/2)^(-a₀ - 2)*gamma(a₀ + 2)
Out[41]:
\begin{equation*}\left(\frac{Y_{1}^{2}}{2} - Y_{1} μ₀ + \frac{Y_{2}^{2}}{2} - Y_{2} μ₀ + \frac{Y_{3}^{2}}{2} - Y_{3} μ₀ + b₀ + \frac{3 μ₀^{2}}{2}\right)^{- a₀ - 2} \Gamma\left(a₀ + 2\right)\end{equation*}
In [42]:
q1_τ = 1/z1_τ * q1_τ
q1_τ
Out[42]:
\begin{equation*}\frac{τ^{a₀ + 1} \left(\frac{Y_{1}^{2}}{2} - Y_{1} μ₀ + \frac{Y_{2}^{2}}{2} - Y_{2} μ₀ + \frac{Y_{3}^{2}}{2} - Y_{3} μ₀ + b₀ + \frac{3 μ₀^{2}}{2}\right)^{a₀ + 2} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ₀ - \frac{Y_{2}^{2}}{2} + Y_{2} μ₀ - \frac{Y_{3}^{2}}{2} + Y_{3} μ₀ - b₀ - \frac{3 μ₀^{2}}{2}\right)}}{\Gamma\left(a₀ + 2\right)}\end{equation*}

数値的に求める

In [43]:
data = [1.1, 1.0, 1.3]
replacements = [(a₀, 1.0), (b₀, 1.0), (μ₀, 0.0), (λ₀, 1.0), zip(Y,data)...]
Out[43]:
7-element Array{Tuple{Sym,Float64},1}:
 (a₀, 1.0)
 (b₀, 1.0)
 (μ₀, 0.0)
 (λ₀, 1.0)
 (Y1, 1.1)
 (Y2, 1.0)
 (Y3, 1.3)
In [44]:
log_p_for_μ_subs = subs(log_p_for_μ, replacements...)
log_p_for_τ_subs = subs(log_p_for_τ, replacements...)
[log_p_for_μ_subs, log_p_for_τ_subs]
Out[44]:
\[ \left[ \begin{array}{r}- 2.0 μ^{2} τ + 3.4 μ τ\\τ \left(- 2.0 μ^{2} + 3.4 μ - 2.95\right) + 2.0 \log{\left (τ \right )}\end{array} \right] \]
In [45]:
q_τ = N(subs(p_τ, replacements...))
q_τ
Out[45]:
\begin{equation*}1.0 e^{- 1.0 τ}\end{equation*}
In [46]:
q_μ = Sym(1)
for i in 1:7
    log_q_μ = N(integrate(log_p_for_μ_subs * q_τ, (τ, 0, oo)))
    z_q_μ = N(integrate(exp(log_q_μ), (μ, -oo, oo)))
    q_μ = 1/z_q_μ * exp(log_q_μ)

    log_q_τ = N(integrate(log_p_for_τ_subs * q_μ, (μ, -oo, oo))(π=>float(π))) # (π=>float(π) が必要なのはちょっと嫌)
    z_q_τ = N(integrate(exp(log_q_τ), (τ, 0, oo)))
    q_τ = 1/z_q_τ * exp(log_q_τ)

    display([q_μ, q_τ])
end
\[ \left[ \begin{array}{r}0.188098154753774 e^{- 2.0 μ^{2} + 3.4 μ}\\4.03007506250001 τ^{2.0} e^{- 2.005 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.112320150163227 e^{- 2.99251870324189 μ^{2} + 5.08728179551122 μ}\\3.11052191637731 τ^{2.0} e^{- 1.83916666666667 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.0965024138432034 e^{- 3.26234707748074 μ^{2} + 5.54599003171727 μ}\\2.97238456804457 τ^{2.0} e^{- 1.81152777777778 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.0938011432750369 e^{- 3.31212144445296 μ^{2} + 5.63060645557004 μ}\\2.94976700750461 τ^{2.0} e^{- 1.8069212962963 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.0933494031271016 e^{- 3.32056521349236 μ^{2} + 5.64496086293701 μ}\\2.94600860516469 τ^{2.0} e^{- 1.80615354938272 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.0932740721319143 e^{- 3.32197669575248 μ^{2} + 5.64736038277921 μ}\\2.94538251532282 τ^{2.0} e^{- 1.80602559156379 τ}\end{array} \right] \]
\[ \left[ \begin{array}{r}0.0932615158350226 e^{- 3.32221205946743 μ^{2} + 5.64776050109463 μ}\\2.94527817564072 τ^{2.0} e^{- 1.80600426526063 τ}\end{array} \right] \]
In [47]:
q_μ*q_τ
Out[47]:
\begin{equation*}0.274681107216063 τ^{2.0} e^{- 1.80600426526063 τ} e^{- 3.32221205946743 μ^{2} + 5.64776050109463 μ}\end{equation*}
In [48]:
julia_code(q_μ*q_τ)
Out[48]:
"0.274681107216063*τ.^2.0.*exp(-1.80600426526063*τ).*exp(-3.32221205946743*μ.^2 + 5.64776050109463*μ)"
In [49]:
"$(q_μ*q_τ)"
Out[49]:
"0.274681107216063*τ^2.0*exp(-1.80600426526063*τ)*exp(-3.32221205946743*μ^2 + 5.64776050109463*μ)"
In [50]:
delta = 0.05
μs = -0.4:delta:2.0
τs = 0.0:delta:5.5

eval(parse("f(μ,τ) = $(q_μ*q_τ)"))
@time c_f = f.(μs', τs);
  0.195528 seconds (718.20 k allocations: 35.829 MiB, 5.98% gc time)
In [51]:
figure(figsize=(5,4))
CS = contour(μs, τs, c_f)
clabel(CS, inline=1, fontsize=10)
xlabel("μ")
ylabel("τ")
grid(ls=":")
In [52]:
figure(figsize=(6.4, 4))
pcolormesh(μs, τs, c_f, cmap="CMRmap")
colorbar()
xlabel("μ")
ylabel("τ")
grid(ls=":")

変分近似無しの計算

分配函数の計算

In [53]:
p3 = simplify(exp(log_p))
Out[53]:
\begin{equation*}\frac{b₀^{a₀} \sqrt{λ₀} τ^{a₀ + 1} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ - \frac{Y_{2}^{2}}{2} + Y_{2} μ - \frac{Y_{3}^{2}}{2} + Y_{3} μ - b₀ - \frac{λ₀ μ^{2}}{2} + λ₀ μ μ₀ - \frac{λ₀ μ₀^{2}}{2} - \frac{3 μ^{2}}{2}\right)}}{4 \pi^{2} \Gamma\left(a₀\right)}\end{equation*}
In [54]:
Z3_μ = integrate(p3, (μ, -oo, oo))
Z3_μ = simplify(Z3_μ)
Out[54]:
\begin{equation*}\frac{\sqrt{2} b₀^{a₀} \sqrt{λ₀} τ^{a₀ + \frac{1}{2}} e^{τ \left(- \frac{Y_{1}^{2}}{2} + \frac{Y_{1}^{2}}{2 \left(λ₀ + 3\right)} + \frac{Y_{1} Y_{2}}{λ₀ + 3} + \frac{Y_{1} Y_{3}}{λ₀ + 3} + \frac{Y_{1} λ₀ μ₀}{λ₀ + 3} - \frac{Y_{2}^{2}}{2} + \frac{Y_{2}^{2}}{2 \left(λ₀ + 3\right)} + \frac{Y_{2} Y_{3}}{λ₀ + 3} + \frac{Y_{2} λ₀ μ₀}{λ₀ + 3} - \frac{Y_{3}^{2}}{2} + \frac{Y_{3}^{2}}{2 \left(λ₀ + 3\right)} + \frac{Y_{3} λ₀ μ₀}{λ₀ + 3} - b₀ + \frac{λ₀^{2} μ₀^{2}}{2 \left(λ₀ + 3\right)} - \frac{λ₀ μ₀^{2}}{2}\right)}}{4 \pi^{\frac{3}{2}} \sqrt{λ₀ + 3} \Gamma\left(a₀\right)}\end{equation*}
In [55]:
F3_μ = expand(-log(Z3_μ))
Out[55]:
\begin{equation*}\frac{Y_{1}^{2} τ}{2} - \frac{Y_{1}^{2} τ}{2 \left(λ₀ + 3\right)} - \frac{Y_{1} Y_{2} τ}{λ₀ + 3} - \frac{Y_{1} Y_{3} τ}{λ₀ + 3} - \frac{Y_{1} λ₀ μ₀ τ}{λ₀ + 3} + \frac{Y_{2}^{2} τ}{2} - \frac{Y_{2}^{2} τ}{2 \left(λ₀ + 3\right)} - \frac{Y_{2} Y_{3} τ}{λ₀ + 3} - \frac{Y_{2} λ₀ μ₀ τ}{λ₀ + 3} + \frac{Y_{3}^{2} τ}{2} - \frac{Y_{3}^{2} τ}{2 \left(λ₀ + 3\right)} - \frac{Y_{3} λ₀ μ₀ τ}{λ₀ + 3} - a₀ \log{\left (b₀ \right )} - a₀ \log{\left (τ \right )} + b₀ τ - \frac{λ₀^{2} μ₀^{2} τ}{2 \left(λ₀ + 3\right)} + \frac{λ₀ μ₀^{2} τ}{2} - \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (τ \right )}}{2} + \frac{\log{\left (λ₀ + 3 \right )}}{2} + \log{\left (\Gamma\left(a₀\right) \right )} + \frac{3 \log{\left (2 \right )}}{2} + \frac{3 \log{\left (\pi \right )}}{2}\end{equation*}
In [56]:
coef = coeff(collect(F3_μ, τ), τ)
Out[56]:
\begin{equation*}\frac{Y_{1}^{2}}{2} - \frac{Y_{1}^{2}}{2 \left(λ₀ + 3\right)} - \frac{Y_{1} Y_{2}}{λ₀ + 3} - \frac{Y_{1} Y_{3}}{λ₀ + 3} - \frac{Y_{1} λ₀ μ₀}{λ₀ + 3} + \frac{Y_{2}^{2}}{2} - \frac{Y_{2}^{2}}{2 \left(λ₀ + 3\right)} - \frac{Y_{2} Y_{3}}{λ₀ + 3} - \frac{Y_{2} λ₀ μ₀}{λ₀ + 3} + \frac{Y_{3}^{2}}{2} - \frac{Y_{3}^{2}}{2 \left(λ₀ + 3\right)} - \frac{Y_{3} λ₀ μ₀}{λ₀ + 3} + b₀ - \frac{λ₀^{2} μ₀^{2}}{2 \left(λ₀ + 3\right)} + \frac{λ₀ μ₀^{2}}{2}\end{equation*}
In [57]:
C = simplify(exp(-simplify(F3_μ - τ*coef))/τ^(a₀+1/Sym(2)))
Out[57]:
\begin{equation*}\frac{\sqrt{2} b₀^{a₀} \sqrt{λ₀}}{4 \pi^{\frac{3}{2}} \sqrt{λ₀ + 3} \Gamma\left(a₀\right)}\end{equation*}
In [58]:
sol1 = solve(diff(coef, Y[1]), Y[1])[1]
Out[58]:
\begin{equation*}\frac{Y_{2} + Y_{3} + λ₀ μ₀}{λ₀ + 2}\end{equation*}
In [59]:
coef1 = simplify(subs(coef, (Y[1], sol1)))
Out[59]:
\begin{equation*}\frac{\frac{Y_{2}^{2} λ₀}{2} + \frac{Y_{2}^{2}}{2} - Y_{2} Y_{3} - Y_{2} λ₀ μ₀ + \frac{Y_{3}^{2} λ₀}{2} + \frac{Y_{3}^{2}}{2} - Y_{3} λ₀ μ₀ + b₀ λ₀ + 2 b₀ + λ₀ μ₀^{2}}{λ₀ + 2}\end{equation*}
In [60]:
sol2 = solve(diff(coef1, Y[2]), Y[2])[1]
Out[60]:
\begin{equation*}\frac{Y_{3} + λ₀ μ₀}{λ₀ + 1}\end{equation*}
In [61]:
coef2 = simplify(subs(coef1, (Y[2], sol2)))
Out[61]:
\begin{equation*}\frac{\frac{Y_{3}^{2} λ₀}{2} - Y_{3} λ₀ μ₀ + b₀ λ₀ + b₀ + \frac{λ₀ μ₀^{2}}{2}}{λ₀ + 1}\end{equation*}
In [62]:
sol = solve(diff(coef2, Y[3]), Y[3])[1]
Out[62]:
\begin{equation*}μ₀\end{equation*}
In [63]:
simplify(subs(coef, ((y, sol) for y in Y)...)) # -> b₀ > 0 and hence coef > 0.
Out[63]:
\begin{equation*}b₀\end{equation*}
In [64]:
ξ = symbols("ξ", positive=true)
Z3 = simplify(integrate(τ^(a₀+1)*exp(-ξ*τ), (τ, 0, oo)))
Z3 = Z3(ξ=>coef)
Z3 = simplify(C*simplify(Z3))
Out[64]:
\begin{equation*}\frac{2^{a₀ + \frac{1}{2}} a₀ b₀^{a₀} \sqrt{λ₀} \left(a₀ + 1\right) \left(λ₀ + 3\right)^{a₀ + \frac{3}{2}} \left(Y_{1}^{2} λ₀ + 2 Y_{1}^{2} - 2 Y_{1} Y_{2} - 2 Y_{1} Y_{3} - 2 Y_{1} λ₀ μ₀ + Y_{2}^{2} λ₀ + 2 Y_{2}^{2} - 2 Y_{2} Y_{3} - 2 Y_{2} λ₀ μ₀ + Y_{3}^{2} λ₀ + 2 Y_{3}^{2} - 2 Y_{3} λ₀ μ₀ + 2 b₀ λ₀ + 6 b₀ + 3 λ₀ μ₀^{2}\right)^{- a₀ - 2}}{\pi^{\frac{3}{2}}}\end{equation*}

事後分布の計算

In [65]:
q3 = p3/Z3
q3 = q3(ξ=>coef) # posterior
Out[65]:
\begin{equation*}\frac{2^{- a₀ - \frac{1}{2}} τ^{a₀ + 1} \left(λ₀ + 3\right)^{- a₀ - \frac{3}{2}} \left(Y_{1}^{2} λ₀ + 2 Y_{1}^{2} - 2 Y_{1} Y_{2} - 2 Y_{1} Y_{3} - 2 Y_{1} λ₀ μ₀ + Y_{2}^{2} λ₀ + 2 Y_{2}^{2} - 2 Y_{2} Y_{3} - 2 Y_{2} λ₀ μ₀ + Y_{3}^{2} λ₀ + 2 Y_{3}^{2} - 2 Y_{3} λ₀ μ₀ + 2 b₀ λ₀ + 6 b₀ + 3 λ₀ μ₀^{2}\right)^{a₀ + 2} e^{τ \left(- \frac{Y_{1}^{2}}{2} + Y_{1} μ - \frac{Y_{2}^{2}}{2} + Y_{2} μ - \frac{Y_{3}^{2}}{2} + Y_{3} μ - b₀ - \frac{λ₀ μ^{2}}{2} + λ₀ μ μ₀ - \frac{λ₀ μ₀^{2}}{2} - \frac{3 μ^{2}}{2}\right)}}{4 \sqrt{\pi} a₀ \left(a₀ + 1\right) \Gamma\left(a₀\right)}\end{equation*}
In [66]:
Nq3 = N(subs(q3, replacements...)) # numerical posterior
Out[66]:
\begin{equation*}\frac{2.41042987827087 τ^{2.0} e^{τ \left(- 2.0 μ^{2} + 3.4 μ - 2.95\right)}}{\sqrt{\pi}}\end{equation*}
In [67]:
simplify(q_μ*q_τ)
Out[67]:
\begin{equation*}0.274681107216063 τ^{2.0} e^{- 3.32221205946743 μ^{2} + 5.64776050109463 μ - 1.80600426526063 τ}\end{equation*}

変分近似との比較するためのプロット

In [68]:
eval(parse("g(μ,τ) = $Nq3")) # numerical posterior
@time c_g = g.(μs', τs);
  0.075039 seconds (225.81 k allocations: 10.485 MiB)
In [69]:
figure(figsize=(10,4))

subplot(121)
CS = contour(μs, τs, c_f)
clabel(CS, inline=1, fontsize=10)
xlabel("μ")
ylabel("τ")
grid(ls=":")
title("variational approximation")

subplot(122)
CS = contour(μs, τs, c_g)
clabel(CS, inline=1, fontsize=10)
xlabel("μ")
ylabel("τ")
grid(ls=":")
title("exact posterior")
Out[69]:
PyObject Text(0.5,1,'exact posterior')
In [70]:
figure(figsize=(10, 3))

subplot(121)
pcolormesh(μs, τs, c_f, cmap="CMRmap")
colorbar()
xlabel("μ")
ylabel("τ")
grid(ls=":")
title("variational approximation")

subplot(122)
pcolormesh(μs, τs, c_g, cmap="CMRmap")
colorbar()
xlabel("μ")
ylabel("τ")
grid(ls=":")
title("exact posterior")
Out[70]:
PyObject Text(0.5,1,'exact posterior')
In [ ]: