階層モデルと予測分布

記事の内容


久しぶりの更新となりました. 今回は階層モデルと予測分布に関して解説します.

問題設定とモデリング

問題設定

次のような問題を考えます. Aさんは15分間に読んだ本のページ数の記録をとっています. ただし, 読む本は観測ごとに異なるものとします. 現在, データが\(N=9\)個ほどあり, 次のような値であるとします.

\[ 11, 10, 9, 8, 14, 13, 12, 13, 11 \]

このデータを用いて, 予測を行いたいと思いました(ここでいう"予測"の意味については, 後ほど). 予測を行うため, 次のようなモデルを仮定しました.

\begin{align} Y_i \mid \lambda_i &\sim \mathrm{Po}(\lambda_i) \quad (i=1,\cdots,N) \\ \lambda_i \mid \beta &\sim \mathrm{Gamma}(\alpha, \beta) \quad (i=1,\cdots, N) \\ \beta &\sim \mathrm{Exp}(\mu_{\beta}) \end{align}

\(Y_i\)は読んだページ数(観測値)とし, \(\alpha=11, \mu_\beta=10^{-4}\)とします. 本ごとに読みやすさ・読みにくさがあることを考慮し, データごとに異なる\(\lambda_i\)を仮定します. また, 次のような記号を用います.

\begin{align} \boldsymbol{Y}_N &= \{ Y_1, \cdots, Y_N\}\\ \boldsymbol{\lambda}_N &= \{ \lambda_1, \cdots, \lambda_N\} \end{align}

ここでの目標は予測分布を計算することです. 解析の前に, "予測"について確認しておきます.

2種類の予測

一般に, 予測分布は次のように定義されます.

【定義: 予測分布】

\(w\in\mathbb{R}^{d}\)をパラメータとする確率モデル\(p(x\mid w)\)に対して, 予測分布を次式で定める. ただし, \(D=\{X_1,\cdots,X_N\}\subset \mathbb{R}^n\)をデータとする.

\begin{equation} p(x\mid D) = \int p(x\mid w)p(w\mid D)dw \end{equation}

今回の問題設定で予測しているものは何でしょうか? 次の2つが考えられます.

  • 各本\(i=1,\cdots,N\)に対して, 新しい観測値\(Y_i^*\)を予測したい.
  • データには含まれない新たな本に対して, 観測値\(Y^*\)を予測したい.

今回はこの2種類の予測のための予測分布を計算します.

解析: 事後分布の計算

本節では, 事後分布を計算します. 事後分布はGibbs samplerにより実現します.

Gibbs samplerの導出

まず, 同時分布を考えます. 観測値\(\boldsymbol{Y}_N\)および, パラメータ\(\boldsymbol{\lambda}_N\), \(\beta\)の同時分布は次のようになります.

\begin{align} p(\boldsymbol{Y}_N, \boldsymbol{\lambda}_N, \beta) &= \left\{ \prod_{i=1}^{N} p(Y_i\mid\lambda_i)p(\lambda_i\mid\beta)) \right\} p(\beta) \\ &\propto \left\{ \prod_{i=1}^{N} \beta^{\alpha}\lambda_i^{y_i+\alpha-1}e^{-(\beta+1)} \right\} e^{-\mu_\beta\beta} \end{align}

したがって, Gibbs samplerで用いる条件付き事後分布は以下のようになります.

\begin{align} p(\lambda_i\mid \boldsymbol{Y}_N, \boldsymbol{\lambda}_N^{\setminus i}, \beta) &\propto p(\boldsymbol{Y}_N, \boldsymbol{\lambda}_N, \beta) \\ &\propto \lambda_i^{y_i+\alpha-1} e^{-(\beta+1)\lambda_i} \\ \\ p(\beta\mid \boldsymbol{Y}_N, \boldsymbol{\lambda}_N) &\propto p(\boldsymbol{Y}_N, \boldsymbol{\lambda}_N, \beta) \\ &\propto \beta^{(N\alpha+1)-1}\mathrm{exp}\left\{ -\left(\mu_{\beta}+\sum_{i=1}^{N}\lambda_i\right) \beta\right\} \end{align}

ここで, \(\boldsymbol{\lambda}_N^{\setminus i}=\boldsymbol{\lambda}_N\setminus \{\lambda_i\}\)とします. 共にガンマ分布です. 次に, この式を用いて事後分布を計算します.

計算とプロット

前のセクションの結果を用いて計算した事後分布からのサンプルのヒストグラムは以下の通りです. 全部で10000サンプルとり, 最初の1000サンプルは初期値に依存するものとして除きました(burn-in期間). 左下の図は\(\boldsymbol{\lambda}_N\)のヒストグラムです. JuliaのPlotsのplot関数で, st=:stephistとすることで, ステップ関数のように描きました. \(\boldsymbol{\lambda}_N\)は観測値の数だけ\(\lambda_i\)を含み, それらを重ねて描いています. 右図は, \(\beta\)の事後分布です.

【コードの実行結果】

次節では, 事後分布からのサンプルを用いて, 予測分布を計算します.

解析: 予測分布の計算

予測分布を計算します. 冒頭でも述べた通り, 階層モデルでは, 予測にも2つの方向性があります. 今回のケースでは, 以下の2つです.

  • 各本\(i=1,\cdots,N\)に対して, 新しい観測値\(Y_i^*\)を予測したい.
  • データには含まれない新たな本に対して, 観測値\(Y^*\)を予測したい.

今回は前者を予測その1, 後者を予測その2として, それぞれ計算します. 予測値であることを強調して*を付けていますが, 以下, 簡単のため\(Y^*\)は単に\(Y\)のように書きます.

予測その1

まずは, 予測その1"各本\(i=1,\cdots,N\)に対して, 新しい観測値\(Y_i\)を予測したい"を考えます. 通常, (教科書的には)予測分布と言えばこちらを指します. この場合, 求める予測分布は以下の通りです.

\begin{equation} p(y_i\mid \boldsymbol{Y}_N) = \int p(y_i\mid \lambda_i)p(\lambda_i\mid \boldsymbol{Y}_N)d\lambda_i \quad (i=1,\cdots, N) \end{equation}

積分が解析的には求まらない場合には, 事後分布からのサンプルを用います. ここでは, 次の2通りの方法で予測分布を近似計算します. 以下, 各\(\lambda_i\)の事後分布からのサンプル\(\{\lambda_i^{(s)}\}_{s=1}^S\)は既に計算してあるものとします.

  • 予測分布の定義式の積分を近似計算し, 確率\(p(y_i\mid \boldsymbol{Y}_N)\)を近似する.
  • 予測分布からのサンプルを乱数発生により得る.

どちらを用いるかは状況次第です. "何が既知か", あるいは"予測分布の何を計算したいか"によります. 前者はモデル内での確率を計算するのに有効です. 一方, 後者は予測分布の平均や分散などを計算するのに有効です.

まずは, 1つ目の計算方法を試します. 確率\(p(y_i\mid \boldsymbol{Y}_N)\)を近似します. 積分は以下のように近似できます.

\begin{align} p(y_i\mid \boldsymbol{Y}_N) &= \int p(y_i\mid \lambda_i)p(\lambda_i\mid \boldsymbol{Y}_N)d\lambda_i \\ &\simeq \frac{1}{S} \sum_{s=1}^{S} p(y_i\mid \lambda_i^{(s)}) \end{align}

ただし, \(\{\lambda_i^{(s)}\}_{s=1}^{S}\)はGibbs samplerによって計算した事後分布からのMCMCサンプルです. これを各\(i=1,\cdots,N\)を固定するごとに, 各\(y_i=0,\cdots,20\)に対して求め, プロットしたものが下図です.

【コードの実行結果】

次に, 予測分布からのサンプルを生成し, 予測分布を構成します. これは以下のようにします. まず, \(i=1,\cdots,N\)を固定します. 次に, 事後分布からのサンプル\(\{\lambda_i^{(s)}\}_{s=1}^{S}\)から, 1つのサンプル\(\lambda_i^{(s)}\)を固定します. そして, モデルにしたがって新たなサンプル \begin{align} y_i^{(s)} &\sim \mathrm{Po}(\lambda_i^{(s)}) \end{align} を計算します. 予測分布は確率モデルを事後分布で重みづけたものですから, このようにして得られたサンプル\(\{y_i^{(s)}\}_{s=1}^{S}\)は予測分布\(p(y_i\mid \boldsymbol{Y}_N)\)からのサンプルと見なすことができます. \(\{y_i^{(s)}\}_{s=1}^{S}\)をヒストグラムで表示したものが下図です.

【コードの実行結果】

予測その2

次に, "データには含まれない新たな本に対して, 観測値\(Y\)を予測したい"を考えます. この場合, 新たな観測値\(Y\)はデータのどの添え字\(i\)にも対応しません. 新たな観測値\(Y\)対応する\(\lambda\)もモデル内にはありません. 元のモデルは以下の通りでした.

\begin{align} Y_i \mid \lambda_i &\sim \mathrm{Po}(\lambda_i) \quad (i=1,\cdots,N) \\ \lambda_i \mid \beta &\sim \mathrm{Gamma}(\alpha, \beta) \quad (i=1,\cdots, N) \\ \beta &\sim \mathrm{Exp}(\mu_{\beta}) \end{align}

これを, 次のようにデータの生成過程とみなします.

\begin{align} Y \mid \lambda &\sim \mathrm{Po}(\lambda) \\ \lambda \mid \beta &\sim \mathrm{Gamma}(\alpha, \beta)\\ \beta &\sim \mathrm{Exp}(\mu_{\beta}) \end{align}

すなわち, \(\beta\)の情報をもとに, 新たに\(\lambda\)を生成し, その\(\lambda\)をもとに観測値\(Y\)を生成します. ただし, 予測分布は事後分布で重みづけられたものですから, 新たな観測値(予測値)は次のように生成されます.

\begin{align} Y \mid \lambda &\sim \mathrm{Po}(\lambda) \\ \lambda \mid \beta &\sim \mathrm{Gamma}(\alpha, \beta)\\ \beta &\sim p(\beta\mid \boldsymbol{Y}_N) \end{align}

\(\beta\)は事後分布から発生させます. この生成過程は次のように解釈できます. まず, 新たな観測値(予測値)\(Y\)に対応する\(\lambda\)を\(\lambda\)の予測分布\(p(\lambda\mid \boldsymbol{Y}_N)\)から生成します. この\(\lambda\)の予測値を用いて, 新たな観測値\(Y\)をモデルにしたがって\(p(y\mid \lambda)\)から生成します. \(Y\)の予測分布を書き下すと, 以下のようになります.

\begin{equation} p(y\mid \boldsymbol{Y}_N) = \int p(y\mid\lambda)p(\lambda\mid \boldsymbol{Y}_N)d\lambda = \int p(y\mid\lambda) \left( \int p(\lambda\mid\beta)p(\beta\mid \boldsymbol{Y}_N) d\beta \right) d\lambda \end{equation}

また, この式を変形すると, 予測分布は次のように書けます.

\begin{align} p(y\mid \boldsymbol{Y}_N) &= \int p(y\mid\lambda) \left( \int p(\lambda\mid\beta)p(\beta\mid \boldsymbol{Y}_N) d\beta \right) d\lambda\\ &= \int \left( \int p(y\mid\lambda)p(\lambda\mid\beta) d\lambda \right) p(\beta\mid \boldsymbol{Y}_N) d\beta \\ \end{align}

すなわち, 予測その2の予測分布は\(\beta\)をパラメータとする確率モデル \begin{equation} p(y\mid\beta) = \int p(y\mid \lambda) p(\lambda\mid \beta)d\lambda \end{equation} に対する予測分布とみなせます(冒頭の予測分布の定義参照).

それでは, 予測分布を計算します. ここでも主に2通りの計算方法が考えられます. 以下, \(\beta\)の事後分布からのサンプル\(\{\beta^{(s)}\}_{s=1}^S\)は既に計算してあるものとします.

  • 予測分布の定義式の積分を近似計算し, 確率\(p(y\mid \boldsymbol{Y}_N)\)を近似する.
  • 予測分布からのサンプルを乱数発生により得る.

まずは1つ目の計算方法からです. 確率\(p(y\mid \boldsymbol{Y}_N)\)を次のように近似します.

\begin{align} p(y\mid \boldsymbol{Y}_N) &= \int p(y\mid\lambda) \left( \int p(\lambda\mid\beta)p(\beta\mid \boldsymbol{Y}_N) d\beta \right) d\lambda\\ &\simeq \frac{1}{S} \sum_{i=1}^{S}p(y\mid \lambda^{(s)})\\ \end{align}

ただし, \(\lambda^{(s)}\)は次のように生成します.

\begin{equation} \lambda^{(s)} \sim \mathrm{Gamma}(\alpha, \beta^{(s)}) \end{equation}

ここで, \(\{\beta^{(s)}\}_{s=1}^{S}\)は事後分布からのサンプルです. このようにして計算した予測分布のプロットが下図です.

【コードの実行結果】

次に, 予測分布からのサンプルを計算して, 予測分布を構成する方法です. 次のように, 生成過程に従って乱数発生させます. すなわち, 各\(s=1,\cdots, S\)に対して, 次のサンプリングを繰り返します.

\begin{align} \lambda^{(s)} &\sim \mathrm{Gamma}(\alpha, \beta^{(s)}) \\ y^{(s)} &\sim \mathrm{Po}(\lambda^{(s)}) \end{align}

ここで, \(\{\beta^{(s)}\}_{s=1}^{S}\)は事後分布からのサンプルです. このようにして構成した予測分布からのサンプル\(\{y^{(s)}\}_{s=1}^{S}\)を下図のようにヒストグラムで表示しました.

【コードの実行結果】

以上, 2通りの計算方法を実行しました. ただし, 今回のケースは, ある程度解析的に処理することもできます. 予測分布を以下のように近似します.

\begin{align} p(y\mid \boldsymbol{Y}_N) &= \int p(y\mid\lambda) \left( \int p(\lambda\mid\beta)p(\beta\mid \boldsymbol{Y}_N) d\beta \right) d\lambda\\ &= \int p(\beta\mid \boldsymbol{Y}_N)\left( \int p(y\mid\lambda)p(\lambda\mid\beta) d\lambda \right)d\beta \\ &= \int p(\beta\mid \boldsymbol{Y}_N) \left\{ \frac{\Gamma (y+\alpha)(\beta^{\alpha}}{y!\Gamma(\alpha)(\beta+1)^{y+\alpha}}\right\} d\beta\\ &\simeq \frac{1}{S}\sum_{i=1}^{S} \frac{\Gamma (y+\alpha)(\beta^{(s)})^{\alpha}}{y!\Gamma(\alpha)(\beta^{(s)}+1)^{y+\alpha}} \end{align}

ここで, \(\{\beta^{(s)}\}_{s=1}^{S}\)は事後分布からのサンプルです. 内部のパラメータに関する積分が解析的に実行できることもあります. 他にも, "数値積分で近似する"などの方法もあります. このようにして近似した予測分布を下図に示します.

【コードの実行結果】

コード

データ等

【Juliaコード1; インポート】
using Plots
using Distributions
using SpecialFunctions
pyplot()
【Juliaコード2; データ】
#データサイズ, サンプル数
n_samps = 10000
Y_data = [11, 10, 9, 8, 14, 13, 12, 13, 11]
N = length(Y_data)

#αとμ_β
α = 11.0
μ_β = 1e-4

#推定
λs, βs = myGibbs_sampler(Y_data, n_samps, α, μ_β)

Gibbs sampler

【Juliaコード3; Gibbs sampler】
function myGibbs_sampler(Y::Array{Int64,1}, n_samps::Int64, α::Float64, μ_β::Float64)
    #burn-in期間
    n_burnin = Int(n_samps/10)

    #サンプルサイズ
    N = length(Y)

    #サンプル保存用
    λsamps = zeros(N, n_samps)
    βsamps = zeros(n_samps)
    λsamps[:,1] = ones(N)
    βsamps[1] = 1.

    #サンプリング
    for s in 2:n_samps
        #λのサンプリング
        for i in 1:N
            λsamps[i,s] = rand(Gamma(Y[i]+α, 1/(βsamps[s-1]+1)))
        end

        #βのサンプリング
        βsamps[s] = rand(Gamma(α*N+1, 1/(μ_β+sum(λsamps[:,s]))))
    end
    return λsamps[:,n_burnin:end], βsamps[n_burnin:end]
end
【Juliaコード4; 事後分布の計算】
#λの事後分布
λ_hist = plot(title="Posterior of λs", xlabel="λ", ylabel="frequency", color=:salmon, legend=false)
for i in 1:N
   plot!(λs[i,:], color=:salmon, st=:stephist, label=false)
end

#βの事後分布
β_hist = plot(βs, title="Posterior of β", xlabel="β", ylabel="frequency", st=:histogram, color=:salmon, legend=false)

#まとめて表示
post_plot = plot(λ_hist, β_hist, size=(1000, 400))
#savefig(post_plot, "figs-HB/post.png")

予測分布

【Juliaコード5; 予測その1】
#モデル
function lik(y::Int64, λ::Float64)
    return pdf(Poisson(λ),y)
end

#予測分布を直接計算
function pred1(y::Int64, λis::Array{Float64,1})
    S = length(λis)
    liks = zeros(S)
    λ = 0.0

    #MC積分
    for s in 1:S
        liks[s] = pdf(Poisson(λis[s]), y)
    end
    return mean(liks)
end

#確率を計算
Y = range(0,20,step=1)
lengthY = length(Y)
preds1 = zeros(N,lengthY)
for i in 1:N
    for j in 1:lengthY
        preds1[i,j] = pred1(Y[j], λs[i,:])
    end
end

#プロット
plot_pred1 = plot(xlabel="Number of pages", ylabel="Probability", title="Predictive distribution", xticks=Y, legend=false)
for i in 1:N
    plot!(Y, preds1[i,:], color=:orange, markershape=:circle, markersize=5, ls=:dash)
end
#savefig(plot_pred1, "figs-HB/pred1.png")
【Juliaコード6; 予測その1】
#予測分布を構成
function pred2(λis::Array{Float64, 1})
    #サンプル数, サンプル保存用
    S = length(λis)
    Ysamps = zeros(S)

    #予測分布からのサンプル
    for s in 1:S
        Ysamps[s] = rand(Poisson(λis[s]))
    end
    return Ysamps
end

#予測分布からのサンプル
preds2 = zeros(N,length(βs))
for i in 1:N
    preds2[i,:] = pred2(λs[i,:])
    println("mean of ith predictive distribution = $(mean(preds2[i,:]))")
end

plot_pred2 = plot(title="Predictive distribution", xlabel="Number of pages", ylabel="frequency", color=:salmon, legend=false, xticks=range(0,20,step=1), xlim=[0,20])
for i in 1:N
   plot!(preds2[i,:], color=:salmon, st=:stephist, label=false, bins=50)
end
#savefig(plot_pred2, "figs-HB/pred2.png")
【Juliaコード7; 予測その2】
#モデル
function lik(y::Int64, λ::Float64)
    return pdf(Poisson(λ), y)
end

#予測分布を直接計算
function pred3(y::Int64, α::Float64, βs::Array{Float64,1})
    S = length(βs)
    liks = zeros(S)
    λ = 0.0

    #MC積分
    for s in 1:S
        #λをサンプル
        λ = rand(Gamma(α, 1/βs[s]))

        #モデル式
        liks[s] = lik(y, λ)
    end
    return mean(liks)
end

#確率を計算
Y = range(0,20,step=1)
lengthY = length(Y)
preds3 = zeros(lengthY)
for j in 1:lengthY
    preds3[j] = pred3(Y[j], α, βs)
end

#プロット
plot_pred3 = plot(xlabel="Number of pages", ylabel="Probability", title="Predictive distribution", xticks=Y, legend=false)
plot!(Y, preds3, color=:orange, markershape=:circle, markersize=5, ls=:dash)
#savefig(plot_pred3, "figs-HB/pred3.png")
【Juliaコー8; 予測その2】
#予測分布を構成
function pred4(α::Float64, βs::Array{Float64,1})
    #サンプル数, サンプル保存用
    S = length(βs)
    Ysamps = zeros(S)
    λ = 0

    #予測分布からのサンプル
    for s in 1:S
        #λをサンプル
        λ = rand(Gamma(α, 1/βs[s]))

        #Yをサンプル
        Ysamps[s] = rand(Poisson(λ))
    end
    return Ysamps
end

#予測分布からのサンプル
preds4 = pred4(α, βs)
println("mean of the predictive distribution = $(round(mean(preds4), digits=3))")
plot_pred4 = plot(preds4, title="Predictive distribution", xlabel="Number of pages", ylabel="probability", st=:histogram, legend=false, normed=true, xticks=range(0,20,step=1), xlim=[0,20])
#savefig(plot_pred4, "figs-HB/pred4.png")
【Juliaコード9; 予測その2】
function lik2(y::Int64, α::Float64, β::Float64)
   return  β^α*gamma(y+α)/gamma(α)/factorial(y)/(β+1)^(y+α)
end

#予測分布を計算
function pred5(y::Int64, α::Float64, βs::Array{Float64,1})
    S = length(βs)
    liks = zeros(S)

    #MC積分
    for s in 1:S
       liks[s] = lik2(y, α, βs[s])
    end
    return mean(liks)
end

#確率を計算
Y = range(0,20,step=1)
lengthY = length(Y)
preds5 = zeros(lengthY)
for j in 1:lengthY
    preds5[j] = pred5(Y[j], α, βs)
end

#プロット
plot_pred5 = plot(xlabel="Number of pages", ylabel="Probability", title="Predictive distribution", xticks=Y, legend=false)
plot!(Y, preds5, color=:orange, markershape=:circle, markersize=5, ls=:dash)
#savefig(plot_pred5, "figs-HB/pred5.png")