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

黒木玄

2018-04-02, 2018-12-20

http://nbviewer.jupyter.org/gist/genkuroki/8a342d0b7b249e279dd8ad6ae283c5db に書いた公式のSymPyによる導出.

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
Y = [symbols("Y$i", real=true) for i in 1:5]
@syms σ λ λ₀ θ α β ξ η positive=true
Out[2]:
(σ, λ, λ₀, θ, α, β, ξ, η)
In [3]:
Y
Out[3]:
\[ \left[ \begin{array}{r}Y_{1}\\Y_{2}\\Y_{3}\\Y_{4}\\Y_{5}\end{array} \right] \]

確率密度函数

In [4]:
p_N = 1/√(2PI)*exp(-λ*(y-μ)^2/2)*√λ
Out[4]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ} e^{- \frac{λ \left(y - μ\right)^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [5]:
p_μ = p_N(μ=>μ₀, y=>μ, λ=>λ*λ₀)
Out[5]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ} \sqrt{λ₀} e^{- \frac{λ λ₀ \left(μ - μ₀\right)^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [6]:
p_λ = exp(-λ/θ)*λ^(α-1)/(gamma(α)*θ^α)
Out[6]:
\begin{equation*}\frac{θ^{- α} λ^{α - 1} e^{- \frac{λ}{θ}}}{\Gamma\left(α\right)}\end{equation*}

対数尤度函数とその二乗の期待値

WAICやWBICを計算するためには対数尤度函数とその二乗の事後分布に関する期待値を計算しなければいけなくなる. 共役事前分布を採用した場合には事後分布も共役事前分布と同じ形になるので、共役事前分布に関する期待値の公式を作っておけば十分である. 以下ではその計算をSymPyを使って行う.

全体を $\mu_0$ だけ平行移動すればいつでも $\mu_0=0$ となるようにできる. 以下ではSymPyでの積分計算を易しくするために $\mu_0=0$ と仮定する.

In [7]:
# for simplicity, put μ₀ = 0.
p_μ = p_μ(μ₀=>0)
Out[7]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ} \sqrt{λ₀} e^{- \frac{λ λ₀ μ^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [8]:
logp_N = simplify(log(p_N))
Out[8]:
\begin{equation*}- \frac{λ \left(y - μ\right)^{2}}{2} + \frac{\log{\left (λ \right )}}{2} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2}\end{equation*}
In [9]:
logp_μ = simplify(log(p_μ))
Out[9]:
\begin{equation*}- \frac{λ λ₀ μ^{2}}{2} + \frac{\log{\left (λ \right )}}{2} + \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2}\end{equation*}
In [10]:
logp_λ = expand(log(p_λ))
Out[10]:
\begin{equation*}- α \log{\left (θ \right )} + α \log{\left (λ \right )} - \log{\left (λ \right )} - \log{\left (\Gamma\left(α\right) \right )} - \frac{λ}{θ}\end{equation*}
In [11]:
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_λ]
Out[11]:
\[ \left[ \begin{array}{r}\frac{\sqrt{2} \sqrt{λ} e^{- \frac{λ \left(y - μ\right)^{2}}{2}}}{2 \sqrt{\pi}}\\\frac{\sqrt{2} \sqrt{λ} \sqrt{λ₀} e^{- \frac{λ λ₀ μ^{2}}{2}}}{2 \sqrt{\pi}}\\\frac{θ^{- α} λ^{α - 1} e^{- \frac{λ}{θ}}}{\Gamma\left(α\right)}\\- \frac{y^{2} λ}{2} + y λ μ - \frac{λ μ^{2}}{2} + \frac{\log{\left (λ \right )}}{2} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2}\\- \frac{λ λ₀ μ^{2}}{2} + \frac{\log{\left (λ \right )}}{2} + \frac{\log{\left (λ₀ \right )}}{2} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2}\\- α \log{\left (θ \right )} + α \log{\left (λ \right )} - \log{\left (λ \right )} - \log{\left (\Gamma\left(α\right) \right )} - \frac{λ}{θ}\end{array} \right] \]
In [12]:
@time I_μ = simplify(integrate(expand(logp_N*p_μ), (μ, -oo, oo)))
  0.667942 seconds (79.46 k allocations: 3.975 MiB)
Out[12]:
\begin{equation*}\frac{λ₀ \left(- y^{2} λ + \log{\left (\frac{λ}{2 \pi} \right )}\right) - 1}{2 λ₀}\end{equation*}
In [13]:
# 対数尤度函数の期待値

@time Elogp_N = simplify(integrate(expand(I_μ*p_λ), (λ, 0, oo)))
  1.858062 seconds (80.13 k allocations: 4.003 MiB)
Out[13]:
\begin{equation*}- \frac{y^{2} α θ}{2} + \frac{\log{\left (θ \right )}}{2} + \frac{\operatorname{polygamma}{\left (0,α \right )}}{2} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2} - \frac{1}{2 λ₀}\end{equation*}
In [14]:
A1 = -2Elogp_N
Out[14]:
\begin{equation*}y^{2} α θ - \log{\left (θ \right )} - \operatorname{polygamma}{\left (0,α \right )} + \log{\left (2 \right )} + \log{\left (\pi \right )} + \frac{1}{λ₀}\end{equation*}
In [15]:
A2 = log(2PI) + α*θ*y^2 + 1/λ₀ - (polygamma(Sym(0),α) + log(θ))
Out[15]:
\begin{equation*}y^{2} α θ - \log{\left (θ \right )} - \operatorname{polygamma}{\left (0,α \right )} + \log{\left (2 \pi \right )} + \frac{1}{λ₀}\end{equation*}
In [16]:
simplify(A1 - A2)
Out[16]:
\begin{equation*}0\end{equation*}
In [17]:
expand(logp_N^2)
Out[17]:
\begin{equation*}\frac{y^{4} λ^{2}}{4} - y^{3} λ^{2} μ + \frac{3 y^{2} λ^{2} μ^{2}}{2} - \frac{y^{2} λ \log{\left (λ \right )}}{2} + \frac{y^{2} λ \log{\left (2 \right )}}{2} + \frac{y^{2} λ \log{\left (\pi \right )}}{2} - y λ^{2} μ^{3} + y λ μ \log{\left (λ \right )} - y λ μ \log{\left (\pi \right )} - y λ μ \log{\left (2 \right )} + \frac{λ^{2} μ^{4}}{4} - \frac{λ μ^{2} \log{\left (λ \right )}}{2} + \frac{λ μ^{2} \log{\left (2 \right )}}{2} + \frac{λ μ^{2} \log{\left (\pi \right )}}{2} + \frac{\log{\left (λ \right )}^{2}}{4} - \frac{\log{\left (\pi \right )} \log{\left (λ \right )}}{2} - \frac{\log{\left (2 \right )} \log{\left (λ \right )}}{2} + \frac{\log{\left (2 \right )}^{2}}{4} + \frac{\log{\left (\pi \right )}^{2}}{4} + \frac{\log{\left (2 \right )} \log{\left (\pi \right )}}{2}\end{equation*}
In [18]:
@time J_μ = integrate(expand(logp_N^2*p_μ), (μ, -oo, oo))
  1.075623 seconds (36 allocations: 896 bytes)
Out[18]:
\begin{equation*}\frac{y^{4} λ^{2}}{4} - \frac{y^{2} λ \log{\left (λ \right )}}{2} + \frac{y^{2} λ \log{\left (2 \right )}}{2} + \frac{y^{2} λ \log{\left (\pi \right )}}{2} + \frac{3 y^{2} λ}{2 λ₀} + \frac{\log{\left (λ \right )}^{2}}{4} - \frac{\log{\left (\pi \right )} \log{\left (λ \right )}}{2} - \frac{\log{\left (2 \right )} \log{\left (λ \right )}}{2} + \frac{\log{\left (2 \right )}^{2}}{4} + \frac{\log{\left (\pi \right )}^{2}}{4} + \frac{\log{\left (2 \right )} \log{\left (\pi \right )}}{2} - \frac{\log{\left (λ \right )}}{2 λ₀} + \frac{\log{\left (2 \right )}}{2 λ₀} + \frac{\log{\left (\pi \right )}}{2 λ₀} + \frac{3}{4 λ₀^{2}}\end{equation*}
In [19]:
# 対数尤度函数の二乗の期待値

@time Elogp_N2 = simplify(integrate(expand(J_μ*p_λ), (λ, 0, oo)))
  8.517231 seconds (24 allocations: 560 bytes)
Out[19]:
\begin{equation*}\frac{y^{4} α^{2} θ^{2}}{4} + \frac{y^{4} α θ^{2}}{4} - \frac{y^{2} α θ \log{\left (θ \right )}}{2} - \frac{y^{2} α θ \operatorname{polygamma}{\left (0,α \right )}}{2} + \frac{y^{2} α θ \log{\left (2 \right )}}{2} + \frac{y^{2} α θ \log{\left (\pi \right )}}{2} + \frac{3 y^{2} α θ}{2 λ₀} - \frac{y^{2} θ}{2} + \frac{\log{\left (θ \right )}^{2}}{4} + \frac{\log{\left (θ \right )} \operatorname{polygamma}{\left (0,α \right )}}{2} - \frac{\log{\left (\pi \right )} \log{\left (θ \right )}}{2} - \frac{\log{\left (2 \right )} \log{\left (θ \right )}}{2} + \frac{\operatorname{polygamma}^{2}{\left (0,α \right )}}{4} - \frac{\log{\left (\pi \right )} \operatorname{polygamma}{\left (0,α \right )}}{2} - \frac{\log{\left (2 \right )} \operatorname{polygamma}{\left (0,α \right )}}{2} + \frac{\operatorname{polygamma}{\left (1,α \right )}}{4} + \frac{\log{\left (2 \right )}^{2}}{4} + \frac{\log{\left (\pi \right )}^{2}}{4} + \frac{\log{\left (2 \right )} \log{\left (\pi \right )}}{2} - \frac{\log{\left (θ \right )}}{2 λ₀} - \frac{\operatorname{polygamma}{\left (0,α \right )}}{2 λ₀} + \frac{\log{\left (2 \right )}}{2 λ₀} + \frac{\log{\left (\pi \right )}}{2 λ₀} + \frac{3}{4 λ₀^{2}}\end{equation*}
In [20]:
@syms ψ₀ ψ₁ 
B1 = 4Elogp_N2(polygamma(Sym(0),α)=>ψ₀, polygamma(Sym(1),α)=>ψ₁)
Out[20]:
\begin{equation*}y^{4} α^{2} θ^{2} + y^{4} α θ^{2} - 2 y^{2} α θ ψ₀ - 2 y^{2} α θ \log{\left (θ \right )} + 2 y^{2} α θ \log{\left (2 \right )} + 2 y^{2} α θ \log{\left (\pi \right )} + \frac{6 y^{2} α θ}{λ₀} - 2 y^{2} θ + ψ₀^{2} + 2 ψ₀ \log{\left (θ \right )} - 2 ψ₀ \log{\left (\pi \right )} - 2 ψ₀ \log{\left (2 \right )} + ψ₁ + \log{\left (θ \right )}^{2} - 2 \log{\left (\pi \right )} \log{\left (θ \right )} - 2 \log{\left (2 \right )} \log{\left (θ \right )} + \log{\left (2 \right )}^{2} + \log{\left (\pi \right )}^{2} + 2 \log{\left (2 \right )} \log{\left (\pi \right )} - \frac{2 ψ₀}{λ₀} - \frac{2 \log{\left (θ \right )}}{λ₀} + \frac{2 \log{\left (2 \right )}}{λ₀} + \frac{2 \log{\left (\pi \right )}}{λ₀} + \frac{3}{λ₀^{2}}\end{equation*}
In [21]:
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(θ))/λ₀)
)
Out[21]:
\begin{equation*}y^{4} α θ^{2} \left(α + 1\right) + \frac{6 y^{2} α θ}{λ₀} - 2 y^{2} \left(α θ \left(ψ₀ + \log{\left (θ \right )}\right) + θ\right) + ψ₁ + \left(ψ₀ + \log{\left (θ \right )}\right)^{2} - 2 \left(ψ₀ + \log{\left (θ \right )}\right) \log{\left (2 \pi \right )} + 2 \left(y^{2} α θ + \frac{1}{λ₀}\right) \log{\left (2 \pi \right )} + \log{\left (2 \pi \right )}^{2} - \frac{2 \left(ψ₀ + \log{\left (θ \right )}\right)}{λ₀} + \frac{3}{λ₀^{2}}\end{equation*}
In [22]:
simplify(B1 - B2)
Out[22]:
\begin{equation*}0\end{equation*}

ベイズ更新

サンプルサイズ $n=5$ の場合の共役事前分布のベイズ公式の公式を求めよう. サンプルによって共役事前分布の4つのパラメーター $\mu_0, \lambda_0, \theta, \alpha$ が更新される.

この節では $\mu_0$ を復活させる.

In [23]:
p_μ = p_μ(μ=>μ-μ₀)
Out[23]:
\begin{equation*}\frac{\sqrt{2} \sqrt{λ} \sqrt{λ₀} e^{- \frac{λ λ₀ \left(μ - μ₀\right)^{2}}{2}}}{2 \sqrt{\pi}}\end{equation*}
In [24]:
n = length(Y)
Ybar = mean(Y)
Out[24]:
\begin{equation*}\frac{Y_{1}}{5} + \frac{Y_{2}}{5} + \frac{Y_{3}}{5} + \frac{Y_{4}}{5} + \frac{Y_{5}}{5}\end{equation*}
In [25]:
Yvar = sum((Y.-Ybar).^2)/length(Y)
Out[25]:
\begin{equation*}\frac{\left(- \frac{Y_{1}}{5} - \frac{Y_{2}}{5} - \frac{Y_{3}}{5} - \frac{Y_{4}}{5} + \frac{4 Y_{5}}{5}\right)^{2}}{5} + \frac{\left(- \frac{Y_{1}}{5} - \frac{Y_{2}}{5} - \frac{Y_{3}}{5} + \frac{4 Y_{4}}{5} - \frac{Y_{5}}{5}\right)^{2}}{5} + \frac{\left(- \frac{Y_{1}}{5} - \frac{Y_{2}}{5} + \frac{4 Y_{3}}{5} - \frac{Y_{4}}{5} - \frac{Y_{5}}{5}\right)^{2}}{5} + \frac{\left(- \frac{Y_{1}}{5} + \frac{4 Y_{2}}{5} - \frac{Y_{3}}{5} - \frac{Y_{4}}{5} - \frac{Y_{5}}{5}\right)^{2}}{5} + \frac{\left(\frac{4 Y_{1}}{5} - \frac{Y_{2}}{5} - \frac{Y_{3}}{5} - \frac{Y_{4}}{5} - \frac{Y_{5}}{5}\right)^{2}}{5}\end{equation*}
In [26]:
# 共役事前分布

prior = simplify(p_μ*p_λ)
Out[26]:
\begin{equation*}\frac{\sqrt{2} θ^{- α} λ^{α - \frac{1}{2}} \sqrt{λ₀} e^{λ \left(- \frac{λ₀ μ^{2}}{2} + λ₀ μ μ₀ - \frac{λ₀ μ₀^{2}}{2} - \frac{1}{θ}\right)}}{2 \sqrt{\pi} \Gamma\left(α\right)}\end{equation*}
In [27]:
p5_N = simplify(prod(p_N(y=>v)^β for v in Y))
Out[27]:
\begin{equation*}\left(\frac{\sqrt{2} λ^{\frac{5}{2}} e^{- \frac{λ \left(\left(Y_{1} - μ\right)^{2} + \left(Y_{2} - μ\right)^{2} + \left(Y_{3} - μ\right)^{2} + \left(Y_{4} - μ\right)^{2} + \left(Y_{5} - μ\right)^{2}\right)}{2}}}{8 \pi^{\frac{5}{2}}}\right)^{β}\end{equation*}
In [28]:
# 次の p5 を分配函数で割ったものが事後分布になる.

p5 = simplify(p5_N*p_μ*p_λ)

# 以下の式を μ, λ に関する定数倍を除いて上の共役事前分布と等しくするためには, 
# 共役事前分布の4つのパラメーターをどのように変えなければいけないかを以下で計算する.
Out[28]:
\begin{equation*}\frac{θ^{- α} λ^{α + \frac{5 β}{2} - \frac{1}{2}} \sqrt{λ₀} \left(2 \pi\right)^{- \frac{5 β}{2} - \frac{1}{2}} e^{λ \left(- \frac{Y_{1}^{2} β}{2} + Y_{1} β μ - \frac{Y_{2}^{2} β}{2} + Y_{2} β μ - \frac{Y_{3}^{2} β}{2} + Y_{3} β μ - \frac{Y_{4}^{2} β}{2} + Y_{4} β μ - \frac{Y_{5}^{2} β}{2} + Y_{5} β μ - \frac{5 β μ^{2}}{2} - \frac{λ₀ μ^{2}}{2} + λ₀ μ μ₀ - \frac{λ₀ μ₀^{2}}{2} - \frac{1}{θ}\right)}}{\Gamma\left(α\right)}\end{equation*}
In [29]:
# n=5 の場合の対数尤度函数
logp5 = expand(log(p5))
Out[29]:
\begin{equation*}- \frac{Y_{1}^{2} β λ}{2} + Y_{1} β λ μ - \frac{Y_{2}^{2} β λ}{2} + Y_{2} β λ μ - \frac{Y_{3}^{2} β λ}{2} + Y_{3} β λ μ - \frac{Y_{4}^{2} β λ}{2} + Y_{4} β λ μ - \frac{Y_{5}^{2} β λ}{2} + Y_{5} β λ μ - α \log{\left (θ \right )} + α \log{\left (λ \right )} - \frac{5 β λ μ^{2}}{2} + \frac{5 β \log{\left (λ \right )}}{2} - \frac{5 β \log{\left (\pi \right )}}{2} - \frac{5 β \log{\left (2 \right )}}{2} - \frac{λ λ₀ μ^{2}}{2} + λ λ₀ μ μ₀ - \frac{λ λ₀ μ₀^{2}}{2} - \frac{\log{\left (λ \right )}}{2} + \frac{\log{\left (λ₀ \right )}}{2} - \log{\left (\Gamma\left(α\right) \right )} - \frac{\log{\left (\pi \right )}}{2} - \frac{\log{\left (2 \right )}}{2} - \frac{λ}{θ}\end{equation*}
In [30]:
# exp() の中の式における -λ の係数
c_λ = -coeff(collect(logp5, λ), λ)
Out[30]:
\begin{equation*}\frac{Y_{1}^{2} β}{2} - Y_{1} β μ + \frac{Y_{2}^{2} β}{2} - Y_{2} β μ + \frac{Y_{3}^{2} β}{2} - Y_{3} β μ + \frac{Y_{4}^{2} β}{2} - Y_{4} β μ + \frac{Y_{5}^{2} β}{2} - Y_{5} β μ + \frac{5 β μ^{2}}{2} + \frac{λ₀ μ^{2}}{2} - λ₀ μ μ₀ + \frac{λ₀ μ₀^{2}}{2} + \frac{1}{θ}\end{equation*}
In [31]:
@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])
$$a = 5 β + λ₀$$
$$b = \frac{Y_{1} β + Y_{2} β + Y_{3} β + Y_{4} β + Y_{5} β + λ₀ μ₀}{5 β + λ₀}$$
$$c = \frac{2 Y_{1}^{2} β^{2} θ + \frac{Y_{1}^{2} β θ λ₀}{2} - Y_{1} Y_{2} β^{2} θ - Y_{1} Y_{3} β^{2} θ - Y_{1} Y_{4} β^{2} θ - Y_{1} Y_{5} β^{2} θ - Y_{1} β θ λ₀ μ₀ + 2 Y_{2}^{2} β^{2} θ + \frac{Y_{2}^{2} β θ λ₀}{2} - Y_{2} Y_{3} β^{2} θ - Y_{2} Y_{4} β^{2} θ - Y_{2} Y_{5} β^{2} θ - Y_{2} β θ λ₀ μ₀ + 2 Y_{3}^{2} β^{2} θ + \frac{Y_{3}^{2} β θ λ₀}{2} - Y_{3} Y_{4} β^{2} θ - Y_{3} Y_{5} β^{2} θ - Y_{3} β θ λ₀ μ₀ + 2 Y_{4}^{2} β^{2} θ + \frac{Y_{4}^{2} β θ λ₀}{2} - Y_{4} Y_{5} β^{2} θ - Y_{4} β θ λ₀ μ₀ + 2 Y_{5}^{2} β^{2} θ + \frac{Y_{5}^{2} β θ λ₀}{2} - Y_{5} β θ λ₀ μ₀ + \frac{5 β θ λ₀ μ₀^{2}}{2} + 5 β + λ₀}{θ \left(5 β + λ₀\right)}$$
In [32]:
# μ に関する平方完成
simplify(c_λ - s[a]/2*(μ-s[b])^2 - s[c])
Out[32]:
\begin{equation*}0\end{equation*}
In [33]:
# α₅ は対数尤度函数における log(λ) の係数に 1/2 を足したもの

α₅ = coeff(collect(logp5, log(λ)), log(λ)) + Sym(1)/2
Out[33]:
\begin{equation*}α + \frac{5 β}{2}\end{equation*}
In [34]:
λ₀₅ = s[a]
Out[34]:
\begin{equation*}5 β + λ₀\end{equation*}
In [35]:
μ₀₅ = s[b]
Out[35]:
\begin{equation*}\frac{Y_{1} β + Y_{2} β + Y_{3} β + Y_{4} β + Y_{5} β + λ₀ μ₀}{5 β + λ₀}\end{equation*}
In [36]:
θ₅ = 1/simplify(s[c])
Out[36]:
\begin{equation*}\frac{θ \left(5 β + λ₀\right)}{2 Y_{1}^{2} β^{2} θ + \frac{Y_{1}^{2} β θ λ₀}{2} - Y_{1} Y_{2} β^{2} θ - Y_{1} Y_{3} β^{2} θ - Y_{1} Y_{4} β^{2} θ - Y_{1} Y_{5} β^{2} θ - Y_{1} β θ λ₀ μ₀ + 2 Y_{2}^{2} β^{2} θ + \frac{Y_{2}^{2} β θ λ₀}{2} - Y_{2} Y_{3} β^{2} θ - Y_{2} Y_{4} β^{2} θ - Y_{2} Y_{5} β^{2} θ - Y_{2} β θ λ₀ μ₀ + 2 Y_{3}^{2} β^{2} θ + \frac{Y_{3}^{2} β θ λ₀}{2} - Y_{3} Y_{4} β^{2} θ - Y_{3} Y_{5} β^{2} θ - Y_{3} β θ λ₀ μ₀ + 2 Y_{4}^{2} β^{2} θ + \frac{Y_{4}^{2} β θ λ₀}{2} - Y_{4} Y_{5} β^{2} θ - Y_{4} β θ λ₀ μ₀ + 2 Y_{5}^{2} β^{2} θ + \frac{Y_{5}^{2} β θ λ₀}{2} - Y_{5} β θ λ₀ μ₀ + \frac{5 β θ λ₀ μ₀^{2}}{2} + 5 β + λ₀}\end{equation*}
In [37]:
α₅ - (α + n*β/2)
Out[37]:
\begin{equation*}0\end{equation*}
In [38]:
λ₀₅ - (λ₀+n*β)
Out[38]:
\begin{equation*}0\end{equation*}
In [39]:
simplify(μ₀₅ - (λ₀*μ₀+n*β*Ybar)/(λ₀+n*β))
Out[39]:
\begin{equation*}0\end{equation*}
In [40]:
tmp = (1/θ)*(1 + (n*β*θ/2)*(λ₀/(λ₀+n*β)*(Ybar - μ₀)^2 + Yvar))
simplify(1/θ₅ - tmp)
Out[40]:
\begin{equation*}0\end{equation*}

まとめ: 以上は $n=5$ の場合の公式. 一般の場合の正規分布モデルの事前共役分布に関するベイズ更新の公式は以下の通り:

$$ \begin{aligned} & \tilde\mu_0 = \frac{\lambda_0\mu_0 + \beta n \bar Y}{\lambda_0 + \beta n}, %\quad & & \tilde\lambda_0 = \lambda_0 + \beta n, \quad \\ & \tilde\alpha = \alpha + \frac{1}{2}\beta n, %\quad & & {\tilde\theta}^{-1} = \theta^{-1}\left[ 1 + \frac{\beta n\theta}{2}\left( \frac{\lambda_0}{\lambda_0 + \beta n}\left(\bar Y - \mu_0\right)^2 + V(Y) \right) \right]. \end{aligned} $$

ここでサイズ $n$ のサンプルを $Y_1,\ldots,Y_n$ と書くと,

$$ \bar Y = \frac{1}{2}\sum_{i=1}^n Y_i, \qquad V(Y) = \frac{1}{2}\sum_{i=1}^n (Y_i - \bar Y)^2. $$

ゆえに $\beta\to\infty$ で

$$ \tilde\mu_0 = \frac{\lambda_0\mu_0 + \beta n \bar Y}{\lambda_0 + \beta n} \to \bar Y, \qquad \tilde\alpha \tilde\theta = \frac{(\alpha + \frac{1}{2}\beta n)\theta} { \left[ 1 + \frac{\beta n\theta}{2}\left( \frac{\lambda_0}{\lambda_0 + \beta n}\left(\bar Y - \mu_0\right)^2 + V(Y) \right) \right] }\to\frac{1}{V(Y)}. $$

このことから $\beta\to\infty$ の極限でベイズ更新は最尤法に一致していることもわかる.

注意: ガンマ分布 $p_\lambda(\lambda|\alpha,\theta)$ の平均と分散はそれぞれ $\alpha\theta$, $\alpha\theta^2$ である. 平均 $\alpha\theta$ を固定したままで分散 $\alpha\theta^2=(\alpha\theta)^2/\alpha$ を小さくするためには $\alpha$ を大きくすればよい. そのベイズ更新結果の平均と分散は $\beta\to\infty$ でそれぞれ $1/V(Y)$ と $0$ に近付く.

正規分布 $p_\mu(\mu|\mu_0,\lambda\lambda_0)$ ($\lambda\lambda_0$ は分散の逆数) の分散を小さくするためには $\lambda\lambda_0$ を大きくすればよい. そのベイズ公式結果の平均と分散はそれぞれ $\bar Y$ と $0$ に近付く.

薄く広がった共役事前分布が欲しければ, $\alpha\theta$ と $\mu_0$ を固定したままで $\alpha$ と $\lambda_0$ を小さくすればよい. 例えば, $\alpha\theta=1$, $\mu_0=0$ と固定したとき, $\alpha$, $\lambda_0$ を小さくして $\theta=1/\alpha$ と設定すれば薄く広がった事前分布が得られる.

事後分布のプロット

In [41]:
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)
Out[41]:
(Y1 => 0.7, Y2 => 1.1, Y3 => 1.2, Y4 => 1.2, Y5 => 1.5)
In [42]:
prior = (p_μ*p_λ)(prior_param...)
Out[42]:
\begin{equation*}\frac{1.58111789863934 \cdot 10^{-8} \sqrt{2} e^{- 1.0 \cdot 10^{-6} λ} e^{- 0.0005 λ μ^{2}}}{\sqrt{\pi} λ^{0.499999}}\end{equation*}
In [43]:
posterior = (p_μ*p_λ)(posterior_param...)(β=>1.0, prior_param..., data_sample...)
posterior = posterior(β=>1.0, prior_param..., data_sample...)
posterior = simplify(posterior)
Out[43]:
\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 [44]:
# 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()
In [45]:
# 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 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 [ ]: