Latent Dirichlet Allocation

記事の内容


久しぶりの新記事. メモも兼ねて, Latent Dirichlet Allocationについて簡単にまとめておきます.

概要

モデルの説明

Latent Dirichlet Allocation(以下,LDA)は, トピックモデルの一種で, 文書のトピック解析に用いられるモデルです. 単語データの生成過程が確率モデルで記述されます. トピックとは, 手元にある文書の"大まかな意味・話題"のようなものです. 例えば手元にあるツイートの集合がおおよそ数学や経済に関するものであると分かれば便利でしょう. また, 手元にあるレビューの集合がある商品に関する否定的意見からなるものであれば, 意思決定に役立つでしょう. 以下, 記号を整理しながらモデルの説明をします.

今, 手元に全部で\(D\)個の文書があるとします. 1つの文書は例えば, 1つのツイートや1つのレビューなどです. 各文書は予め単語で分けてあり, 全\(V\)種類あるとします. 単語には\(1\)から\(V\)までの番号をつけておき, 以降, 単語とこの番号を同一視します. \(d\)番目の文書の\(i\)番目の単語を\(w_{di}\)で表します. また, \(d\)番目の文書の単語数を\(N_d\)で表します.

一般に, 単語には出現しやすい特定の文脈があります. 例えば, "確率"や"統計"などの単語は数学や工学などのトピックに現れやすく, 芸能やスポーツなどのトピックには現れにくいはずです. そこで, 各単語に対してどのトピックに属するのかを表現する変数を導入します. トピック数\(K\)は予め固定した上で, 単語\(w_{di}\)に対するトピックを表す変数を\(z_{di}\)で表します. \(z_{di}\)は\(1\)から\(K\)までの整数値をとります. 以下, トピックとトピックの番号を同一視します.

次に, 単語とトピックの出やすさを確率で表現します. 単語の出やすさはトピックに大きく依存します. そこで, \(k\)番目のトピックにおける単語の出現確率をベクトル\(\boldsymbol{\phi}_k\)に格納します. 一方, トピックの出やすさは文書に大きく依存します. そこで\(d\)番目の文書におけるトピック\(k\)の出現確率をベクトル\(\boldsymbol{\theta}_d\)に格納します.

以上の考察を踏まえて, 確率分布を用いてデータの生成過程を次のようにモデル化します. このモデルでは, \(\boldsymbol{\alpha}\), \(\boldsymbol{\beta}\)とトピック数\(K\)を予め与えます. また, 文書数\(D\), 文書長\(N_d\), 単語種類数\(V\), 単語はデータから与えられます.

  • 各トピック\(k=1,\cdots,K\)に対して, 単語の出現確率\(\boldsymbol{\phi}_k\in\mathbb{R}^V\)がDirichlet分布に従って得られる.
  • \begin{equation} \boldsymbol{\phi}_k \sim \mathrm{Dir}(\boldsymbol{\beta}) \end{equation}
  • 各文書\(d=1,\cdots,D\)に対して, トピックの出現確率\(\boldsymbol{\theta}_d\in\mathbb{R}^K\)がDirichlet分布に従って得られる.
  • \begin{equation} \boldsymbol{\theta}_d \sim \mathrm{Dir}(\boldsymbol{\alpha}) \end{equation}
  • 各文書\(d=1,\cdots,D\)の\(i=1,\cdots,N_d\)番目の単語に対応するトピック\(z_{di}\)が確率分布\(\boldsymbol{\theta}_d\)に従って得られる.
  • \begin{equation} p(z_{di}\mid\boldsymbol{\theta}_d) = \theta_{z_{di},d} \end{equation}
  • 各文書\(d=1,\cdots,D\)の\(i=1,\cdots,N_d\)番目の単語\(w_{di}\)が確率分布\(\boldsymbol{\phi}_{z_{di}}\)に従って得られる.
  • \begin{equation} p(w_{di}\mid\boldsymbol{\phi}_{z_{di}}, z_{di}) = \phi_{w_{di}, z_{di}} \end{equation}

記号をいくつか導入しておきます. 単語とトピックを次のようにまとめておきます.

\begin{equation} \boldsymbol{w}_d = \{w_{di}\}_{i=1}^{N_d},\quad \boldsymbol{z}_d=\{z_{di}\}_{i=1}^{N_d} \end{equation}

また, 確率変数を以下のようにまとめておきます.

\begin{equation} W=\{\boldsymbol{w}_1, \cdots,\boldsymbol{w}_D\},\quad Z=\{\boldsymbol{z}_1, \cdots,\boldsymbol{z}_D\},\quad \Theta=\{\boldsymbol{\theta}_1, \cdots,\boldsymbol{\theta}_D\},\quad \Phi=\{\boldsymbol{\phi}_1, \cdots,\boldsymbol{\phi}_D\} \end{equation}

LDAをモデルとした推論では, \(W\)から\(Z,\Theta,\Phi\)を推論するのが目標です.

ちなみに同時分布を書き下すと以下のようになります. 複雑な形なので近似推論法を使うのが自然です.

\begin{align} p(W,Z,\Theta,\Phi) &= \left\{ \prod_{d=1}^{D}\prod_{i=1}^{N_d} p(w_{di}\mid z_{di},\boldsymbol{\phi}_{z_{di}})p(z_{di}\mid\boldsymbol{\theta}_d)\right\}\left\{\prod_{d=1}^D p(\boldsymbol{\theta}_d)\right\}\left\{\prod_{k=1}^K p(\boldsymbol{\phi}_k )\right\} \\ &\propto \left\{ \prod_{d=1}^{D}\prod_{i=1}^{N_d} \phi_{w_{di},z_{di}}\theta_{z_{di},d}\right\}\left\{\prod_{d=1}^D \left(\prod_{k=1}^K\theta_{kd}^{\alpha_k-1}\right)\right\}\left\{\prod_{k=1}^K \left(\prod_{v=1}^V\phi_{vk}^{\beta_v-1}\right)\right\} \end{align}

データの生成過程

先ほどのモデルをアルゴリズムっぽく書いておきます. 実際の分析では, トピック数\(K\), 事前分布のパラメータ\(\boldsymbol{\alpha},\boldsymbol{\beta}\)は自分で設定します. 一方, 単語データおよび, 文書数\(D\), 文書長\(N_d\), 単語種類数\(V\)はデータから分かります. また, 具体的なトピックやトピックの順番は指定せず, トピック数のみを指定します.

データの生成過程は次のように記述できます.

【LDAのデータ生成】

Set the global parameter \(\boldsymbol{\alpha}\), \(\boldsymbol{\beta}\) and the number of topics \(K\).
for \(k=1,\cdots,K\)
Sample \(\boldsymbol{\phi}_k \sim \mathrm{Dir}(\boldsymbol{\beta})\).
end
for \(d=1,\cdots,D\)
Sample \(\boldsymbol{\theta}_d \sim \mathrm{Dir}(\boldsymbol{\alpha})\).
for \(i=1,\cdots,N_d\)
Sample \(z_{di}\) with probability \(p(z_{di}\mid\boldsymbol{\theta}_d)=\theta_{z_{di},d}\)
Sample \(w_{di}\) with probability \(p(w_{di}\mid\boldsymbol{\phi}_{z_{di}})=\phi_{w_{di}, z_{di}}\)
end
end

【実験】データを作ってみる

モデルから実際に人工的にデータを生成してみます. 文章数\(D=100\), 単語種類数\(V=700\), トピック数\(K=10\)とし, 文書長\(N_d(d=1,\cdots,D)\)をPoisson分布\(\mathrm{Poisson}(10)\)から発生させます. また, \(\boldsymbol{\alpha}=\frac{2}{K}\mathrm{rand}(K)\), \(\boldsymbol{\beta}=\frac{10}{V}\mathrm{rand}(V)\)とします.

各文書でのトピック出現確率\(\boldsymbol{\theta}_d\)の成分を棒グラフで示します. \(\boldsymbol{\theta}_d\)は全部で\(D=100\)個あります. 下図に\(d=10,20,\cdots,100\)のみを抜き出しました. 例えば, \(10\)番目の文書には\(9\)番目のトピックが現れやすいことが分かります.

【コード4の実行結果】

下図は, 各トピックの単語出現確率を棒グラフで示したものです.

【コード5の実行結果】

以下の2節では, このデータを用いて推論を行います. 推論方法として, Gibbs samplerと変分推論を用います.

推論1: Gibbs sampler

更新式の導出

まずはGibbs samplerを用いた推論です. 条件付き事後分布を求めると次のようになります.

\begin{align} z_{di}\mid W,Z^{\setminus {di}}, \Theta, \Phi &\sim \mathrm{Categorical}(\boldsymbol{\hat{\boldsymbol{p}}}^{(di)}) ,\quad \hat{p}^{(di)}_k \propto \phi_{w_{di},k}\theta_{kd} \\ \boldsymbol{\theta}_d \mid W, Z, \Theta^{\setminus d}, \Phi &\sim \mathrm{Dir}(\hat{\boldsymbol{\alpha}}^{(d)}) ,\quad \hat{\alpha}^{(d)}_k = \alpha_k + \sum_{i=1}^{N_d}\mathbb{I}[z_{di}=k] \\ \boldsymbol{\phi}_k \mid W, Z, \Theta, \Phi^{\setminus k} &\sim \mathrm{Dir}(\hat{\boldsymbol{\beta}}^{(k)}), \quad \hat{\beta}^{(k)}_v = \beta_v + \sum_{v=1}^V\sum_{i=1}^{N_d}\mathbb{I}[w_{di}=v]\mathbb{I}[z_{di}=k] \end{align}

アルゴリズムにまとめると以下のようになります.

【LDAに対するGibbs sampler】

Initialize \(Z^{(0)}\), \(\Theta^{(0)}\), \(\Phi^{(0)}\)
for \(s=1,2,\cdots\)
for \(d=1,\cdots,D\)
for \(i=1,\cdots,N_d\)
Sample \(z_{di}^{(s)}\) from \(\mathrm{Categorical}(\hat{\boldsymbol{p}}^{(di)})\)
end
Sample \(\boldsymbol{\theta}_d^{(s)} \) from \(\mathrm{Dir}(\hat{\boldsymbol{\alpha}}^{(d)})\)
end
for \(k=1,\cdots,K\)
Sample \(\boldsymbol{\phi}_k^{(s)}\) from \(\mathrm{Dir}(\hat{\boldsymbol{\beta}}^{(k)})\)
end
end

【実験】推論してみる

上のアルゴリズムを用いて実験します. データは前節で生成したものを用い, モデルもLDAです. データから文書数\(D\)や文書長\(N_d\), 単語の種類数\(V\)は分かります. 自分で設定するのはトピック数\(K\)とパラメータ\(\boldsymbol{\alpha}\), \(\boldsymbol{\beta}\)です. トピック数は\(K=10\)とし, \(\boldsymbol{\alpha}=\frac{1}{K}\mathrm{ones}(K)\), \(\frac{100}{V}\mathrm{ones}(V)\)としました. サンプル数は\(10000\)とし, そのうち\(1000\)サンプルを初期値に依存するものとして除去しました.

下図に\(\boldsymbol{\theta}_d\)の値を棒グラフで示します. トピック番号の並びは関係ないので, 真値とは棒グラフの様子も異なります.

【コード8の実行結果】

下図に\(\boldsymbol{\phi}_k\)の値を棒グラフで示します.

【コード9の実行結果】

推論2: 変分推論

更新式の導出

変分推論を用います. ここでは平均場近似を用います. 次の分布で事後分布を近似します.

\begin{equation} r(Z,\Theta,\Phi) = \left( \prod_{d=1}^D\prod_{i=1}^{N_d}r(z_{di})\right)\left( \prod_{d=1}^Dr(\boldsymbol{\theta}_d)\right) \left( \prod_{k=1}^Kr(\boldsymbol{\phi}_k)\right) \end{equation}

平均場近似の近似分布の式から, 次式が得られます.

\begin{align} \boldsymbol{\theta}_d &\sim \mathrm{Dir}(\hat{\boldsymbol{\alpha}}^{(d)}),\quad \hat{\alpha}^{(d)}_k = \alpha_k + \sum_{i=1}^{N_d}r(z_{di}=k) \\ \boldsymbol{\phi}_k &\sim \mathrm{Dir}(\hat{\boldsymbol{\beta}}^{(k)}),\quad \hat{\beta}^{(k)}_v = \beta_v + \sum_{d=1}^D\sum_{i=1}^{N_d}\mathbb{I}[w_{di}=v]r(z_{di}=k) \\ z_{di} &\sim \mathrm{Categorical}(\hat{\boldsymbol{p}}^{(di)}),\quad \hat{p}^{(di)}_k \propto \exp\left\{ \Psi\left(\hat{\alpha}_{z_{di}}^{(d)}\right)+\Psi\left(\hat{\beta}_{w_{di}}^{(k)}\right)-\Psi\left(\sum_{k=1}^K\hat{\alpha}_k^{(d)}\right)-\Psi\left(\sum_{v=1}^V\hat{\beta}_v^{(k)}\right)\right\} \end{align}

以下にアルゴリズムを示します. 実装方法は様々あると思います. 以下のアルゴリズムでは, \(\alpha^{(d)}_k=\alpha_k+\Delta_{kd}\)のように, 更新により足される部分を\(\Delta\)を用いて表しています. 反復ごとにこれを計算しなおします. ただし最初は乱数でランダムに初期化します.

【LDAに対する変分推論】

Initialize \(\Delta^{(kd)}\), \(\Delta^{(vk)}\) randomly and Initialize \(\boldsymbol{p}^{(di)}\).
for \(n=1,2,\cdots\)
for \(d=1,\cdots,D\)
for \(i=1,\cdots,N_d\)
for \(k=1,\cdots,K\)
\(\Delta^{(kd)}+=r(z_{di}=k)\)
\(\Delta^{(w_{di}k)}+=r(z_{di}=k)\)
Update \(p^{(di)}_k\) using \(\hat{\boldsymbol{\alpha}}^{(d)}\), \(\hat{\boldsymbol{\beta}}^{(k)}\).
end
end
end
Update \(\hat{\alpha}^{(d)}_k=\alpha_k+\Delta^{(kd)}\), \(\hat{\beta}^{(k)}_v = \beta_v+\Delta^{(kv)}\).
Initialize \(\Delta^{(kd)}=0, \Delta^{(vk)}=0\).
end

【実験】推論してみる

Gibbs samplerの場合と設定は同じです. 反復回数を\(1000\)としました. 以下に推定した\(\boldsymbol{\theta}_d\)を棒グラフで示します. こちらもトピックの順番は関係ないのでまたもや様子は異なります.

【コード12の実行結果】

下図は\(\boldsymbol{\phi}_k\)の推定結果です.

【コード13の実行結果】

コード

【Juliaコード1; 初期化】
#mathematics
using LinearAlgebra
using SpecialFunctions

#statistics
using Random
using Statistics
using Distributions

#visualize
using Plots
pyplot()

#macros
using UnPack
using ProgressMeter
【Juliaコード2; データ生成】
function create_data(αvec, βvec, K, V, D, Nds)
    n_words = sum(Nds)
    Θ = zeros(K, D) #topic probability θd vectors
    Z = zeros(Int64, n_words) #topic vector
    Φ = zeros(V, K) #topic probability ϕk vectors
    W = zeros(Int64, n_words) #word vector
    
    #create topic 
    for k in 1:K
        Φ[:,k] = rand(Dirichlet(βvec))
    end
    
    #create words
    idx = 0
    for d in 1:D
        Θ[:,d] = rand(Dirichlet(αvec))
        for i in 1:Nds[d]
            idx += 1
            Z[idx] = rand(Categorical(Θ[:,d]))
            W[idx] = rand(Categorical(Φ[:,Z[idx]]))
        end
    end
    return W, Z, Θ, Φ
end
【Juliaコード3; データ生成】
#set the random seed
Random.seed!(42)

#paramters
D = 100
K = 10
V = 700
Nds = rand(Poisson(10), D)
αvec = 2*rand(K)/K
βvec = 10*rand(V)/V

#data
W, Z_true, Θ_true, Φ_true = create_data(αvec, βvec, K, V, D, Nds)
data = (W=W, D=D, V=V, Nds=Nds)
【Juliaコード4; トピック出現確率】
#θ
title = plot(title="θ:topic probability", grid=false, showaxis = false, ticks = false)
local_titles = ["θ₁₀", "θ₂₀", "θ₃₀", "θ₄₀", "θ₅₀", "θ₆₀", "θ₇₀", "θ₈₀", "θ₉₀", "θ₁₀₀"]
plot_array = []
for i in 1:10
    p = plot(1:K, Θ_true[:,10*i], st=:bar, label=false, ylim=[0,1], xticks=1:K, title=local_titles[i]
        , xlabel="topic", ylabel="probability")
    push!(plot_array, p)
end
fig1 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig1.png")
【Juliaコード5; 単語出現確率】
#ϕ
title = plot(title="ϕ:word probability", grid=false, showaxis = false, ticks = false)
local_titles = ["ϕ₁", "ϕ₂", "ϕ₃", "ϕ₄", "ϕ₅", "ϕ₆", "ϕ₇", "ϕ₈", "ϕ₉", "ϕ₁₀"]
plot_array = []
for i in 1:10
    p = plot(1:V, Φ_true[:,i], st=:bar, label=false, ylim=[0,1], xticks=0:200:V, title=local_titles[i]
    , xlabel="word", ylabel="probability")
    push!(plot_array, p)
end
fig2 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig2.png")
【Juliaコード6; Gibbs sampler】
#initialize the probability vectors
function init_probs(αvec, βvec, D, V, K, n_samps)
    Θsamps = zeros(K, D, n_samps)
    Φsamps = zeros(V, K, n_samps) 
    for d in 1:D
        Θsamps[:,d,1] = rand(Dirichlet(αvec))
    end
    for k in 1:K
        Φsamps[:,k,1] = rand(Dirichlet(βvec))
    end
    return Θsamps, Φsamps
end

#update probability vector p from p1 and p2
function update_pvec(p1, p2)
    tmp = p1 .* p2
    return tmp/sum(tmp)
end

#update counter
function count(WZcount, Zcount, d, w, z)
    WZcount[w,z] += 1
    Zcount[z,d] += 1
end

function my_Gibbs_sampler(data, model_params, n_samps, n_burnin)
    @unpack W,D,V,Nds = data
    @unpack K,αvec,βvec = model_params
    
    #initialize
    n_words = sum(Nds)
    Θsamps, Φsamps = init_probs(αvec, βvec, D, V, K, n_samps)
    Zsamps = zeros(Int64, n_words, n_samps)
    pvec = zeros(K)
    
    #counter
    Zcount = zeros(K,D)
    WZcount = zeros(V,K)
    
    #sample loop
    idx = 0
    @inbounds @showprogress for s in 2:n_samps
        for d in 1:D
            for i in 1:Nds[d]
                idx += 1
                pvec = update_pvec(Φsamps[W[idx],:,s-1], Θsamps[:,d,s-1])
                Zsamps[idx,s] = rand(Categorical(pvec))
                count(WZcount, Zcount, d, W[idx], Zsamps[idx,s])
            end
            Θsamps[:,d,s] = rand(Dirichlet(αvec + Zcount[:,d]))
        end
        for k in 1:K
            Φsamps[:,k,s] = rand(Dirichlet(βvec + WZcount[:,k]))
        end
        idx = 0
    end
    return Θsamps[:,:,n_burnin:end], Φsamps[:,:,n_burnin:end], Zsamps[:,n_burnin:end]
end
【Juliaコード7; Gibbs sampler実行】
#model paramters
Random.seed!(42)
K = 10
αvec = ones(K)/K
βvec = 100*ones(V)/V
model_params = (K=K, αvec=αvec, βvec=βvec)

#Gibbs sampler
n_samps = 10000
n_burnin = div(n_samps, 10)
@time Θsamps, Φsamps, Zsamps = my_Gibbs_sampler(data, model_params, n_samps, n_burnin)

#estimates
Θest = mean(Θsamps, dims=3)
Φest = mean(Φsamps, dims=3)
Zest = Zsamps[:,end]
【Juliaコード8; トピック出現確率】
#θ
title = plot(title="θ:topic probability", grid=false, showaxis = false, ticks = false)
local_titles = ["θ₁₀", "θ₂₀", "θ₃₀", "θ₄₀", "θ₅₀", "θ₆₀", "θ₇₀", "θ₈₀", "θ₉₀", "θ₁₀₀"]
plot_array = []
for i in 1:10
    p = plot(1:K, Θest[:,10*i], st=:bar, label=false, ylim=[0,1], xticks=1:K, title=local_titles[i]
    , xlabel="topic", ylabel="probability")
    push!(plot_array, p)
end
fig3 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig3.png")
【Juliaコード9; 単語出現確率】
#ϕ
title = plot(title="ϕ:word probability", grid=false, showaxis = false, ticks = false)
local_titles = ["ϕ₁", "ϕ₂", "ϕ₃", "ϕ₄", "ϕ₅", "ϕ₆", "ϕ₇", "ϕ₈", "ϕ₉", "ϕ₁₀"]
plot_array = []
for i in 1:10
    p = plot(1:V, Φest[:,i], st=:bar, label=false, ylim=[0,1], xticks=0:200:V, title=local_titles[i]
    , xlabel="word", ylabel="probability")
    push!(plot_array, p)
end
fig4 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig4.png")
【Juliaコード10; 変分推論】
#Dirichlet expectation
ElogDir(vec, i) = digamma(vec[i])-digamma(sum(vec))
    
ELBO(pvechats, αvechats, βvechats, αvec, βvec, W, Nds, K, D, ELBO_const) = (
    ELBO_const - sumlogconstDir(αvechats, βvechats, D, K)
    + Evarloglik(pvechats, αvechats, βvechats, W, Nds, D, K)
    + Evarlogpθ(αvechats, αvec, D, K)
    + Evarlogpϕ(βvechats, βvec, V, K)
    - Ezlogrz(pvechats, Nds, D, K)
    - EvarlogrΘ(αvechats, D, K)
    - EvarlogrΦ(βvechats, V, K)
)   

#update probability
function update_prob(αvechats, βvechats, w, d, K)
    pvec = zeros(K)
    for k in 1:K
        pvec[k] = exp(ElogDir(αvechats[:,d], k) + ElogDir(βvechats[:,k], w)) 
    end
    return pvec/sum(pvec)
end

#variational infernce
function myVI(data, model_params, n_train)
    @unpack W,D,V,Nds = data
    @unpack K,αvec,βvec = model_params
    #initialize
    w = 0
    idx = 0
    n_words = sum(Nds)
    αvecs = αvec*ones(D)'
    βvecs = βvec*ones(K)'
    Eθvecs = rand(K,D)
    Eϕvecs = rand(V,K)
    αvechats = αvecs + Eθvecs
    βvechats = βvecs + Eϕvecs
    pvechats = ones(K)/K * ones(n_words)'
    @showprogress for n in 1:n_train
        idx = 0
        for d in 1:D
            for i in 1:Nds[d]
                idx += 1
                w = W[idx]
                for k in 1:K
                    Eθvecs[k,d] += pvechats[k,idx]
                    Eϕvecs[w,k] += pvechats[k,idx]
                end
                pvechats[:,idx] = update_prob(αvechats, βvechats, w, d, K)
            end
        end
        αvechats = αvecs + Eθvecs
        βvechats = βvecs + Eϕvecs
        Eθvecs = zeros(K,D)
        Eϕvecs = zeros(V,K)
    end
    return pvechats, αvechats, βvechats
end

#estimate
function estimate(paramvecs)
    dim,N = size(paramvecs)
    Eparams = zeros(dim,N)
    for n in 1:N
        Eparams[:,n] = paramvecs[:,n]/sum(paramvecs[:,n])
    end
    return Eparams
end
【Juliaコード11; 変分推論の実行】
#model paramters
Random.seed!(42)
K = 10
αvec = ones(K)/K
βvec = 100*ones(V)/V
model_params = (K=K, αvec=αvec, βvec=βvec)

#variational inference
n_train = 1000
@time pvechats, αvechats, βvechats = myVI(data, model_params, n_train)

#estimates(mean)
Θest = estimate(αvechats)
Φest = estimate(βvechats)
【Juliaコード12; トピック出現確率】
#θ
title = plot(title="θ:topic probability", grid=false, showaxis = false, ticks = false)
local_titles = ["θ₁₀", "θ₂₀", "θ₃₀", "θ₄₀", "θ₅₀", "θ₆₀", "θ₇₀", "θ₈₀", "θ₉₀", "θ₁₀₀"]
plot_array = []
for i in 1:10
    p = plot(1:K, Θest[:,10*i], st=:bar, label=false, ylim=[0,1], xticks=1:K, title=local_titles[i]
    , xlabel="topic", ylabel="probability")
    push!(plot_array, p)
end
fig5 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig5.png")
【Juliaコード13; 単語出現確率】
#ϕ
title = plot(title="ϕ:word probability", grid=false, showaxis = false, ticks = false)
local_titles = ["ϕ₁", "ϕ₂", "ϕ₃", "ϕ₄", "ϕ₅", "ϕ₆", "ϕ₇", "ϕ₈", "ϕ₉", "ϕ₁₀"]
plot_array = []
for i in 1:10
    p = plot(1:V, Φest[:,i], st=:bar, label=false, ylim=[0,1], xticks=0:200:V, title=local_titles[i]
    , xlabel="word", ylabel="probability")
    push!(plot_array, p)
end
fig6 = plot(title, plot_array..., layout=@layout[a{0.01h}; grid(2,5)], size=(900,400))
savefig("figs-LDA/fig6.png")
参考文献

      [1]佐藤一誠, トピックモデルによる統計的潜在意味解析, 2018, コロナ社
      [2]岩田具治, トピックモデル, 2015, 講談社