都内の某IT企業にてデータサイエンティストとして働いている @yusaku-i です.
メイン言語を から に移行して半年ほど経ちます.
ローカル環境で粛々と開発してきた機械学習手法たちを整理してGitHubにパッケージとして公開しています.

今回は確率的ブロックモデル(SBM; stochastic block model)をスクラッチ実装したので紹介したいと思います.
体系的に学びたい方は「関係データ学習(機械学習プロフェッショナルシリーズ)」を参考にしてください.

SBMとは

関係データはSNSでの「フォロー・友達」やウェブサイトでの「ハイパーリンク」,ECサイトでの「購買履歴」など,2つのオブジェクト間にリンクが存在するか否かで定義されるデータになります.
なお,本記事では「関係の有無」は2値の離散データのみで表されるとします(連続値を用いると関係の強さなどを考慮できるみたいです).

SBMは関係データ行列が潜在的なブロック構造をもつことを仮定した確率的生成モデルです.

286eb0b0-a745-438c-b9f0-532837c7a3dd.png

SBMを使うことで右図のように行と列をクラスタリングして「リンクが存在しやすい」ブロックと「リンクが存在しにくい」ブロックとで表現することができます.

定式化

関係データ行列の行 $i=1,\ldots,N^\textit{1}$ を第1ドメイン,列 $j=1,\ldots,N^\textit{2}$ を第2ドメインとし,それぞれ $K$ 個と $L$ 個の潜在クラスタをもつとします.
このときSBMは次のような生成過程をもちます.

\begin{align*}
{\boldsymbol \theta}^\textit{1} &\sim {\rm Dirichlet}({\boldsymbol \alpha}^\textit{1}) \\
{\boldsymbol \theta}^\textit{2} &\sim {\rm Dirichlet}({\boldsymbol \alpha}^\textit{2}) \\
\{\{ \phi_{k\ell} \}_{\ell=1}^{L} \}_{k=1}^{K} &\sim {\rm Beta}(a_0, b_0) \\
\{ z_i^\textit{1} \}_{i=1}^{N^\textit{1}} &\sim {\rm Categorical}({\boldsymbol \theta}^\textit{1}) \\
\{ z_j^\textit{2} \}_{j=1}^{N^\textit{2}} &\sim {\rm Categorical}({\boldsymbol \theta}^\textit{2}) \\
\{\{ x_{ij} \}_{j=1}^{N^\textit{2}} \}_{i=1}^{N^\textit{1}} &\sim {\rm Bernoulli}(\phi_{z_i^\textit{1} z_j^\textit{2}})
\end{align*}

${\boldsymbol \theta}$ はクラスタの混合割合を意味し,行と列とでそれぞれディリクレ分布を事前分布とします.
また,${\boldsymbol \Phi}=\{\phi _{k\ell}\}$ はリンクが存在するかの $x \in \{0,1 \}$ を表すベルヌーイ分布であり,ベータ分布を事前分布とします.
ここで,ディリクレ分布やベータ分布は共役事前分布と呼ばれる確率分布です.

潜在クラスタを ${\bf z}$ とするとSBMの同時確率は次のようになります.

\begin{align*}
&p({\bf x}, {\bf z}^\textit{1}, {\bf z}^\textit{2}, {\boldsymbol \theta}^\textit{1}, {\boldsymbol \theta}^\textit{2}, {\boldsymbol \Phi}\ \vert \ {\boldsymbol \alpha}^\textit{1}, {\boldsymbol \alpha}^\textit{2}, a_0, b_0) \\
&= p({\bf z}^\textit{1}\ \vert \ {\boldsymbol \theta}^\textit{1}) \cdot p({\boldsymbol \theta}^\textit{1}\ \vert \ {\boldsymbol \alpha}^\textit{1}) \cdot p({\bf z}^\textit{2}\ \vert \ {\boldsymbol \theta}^\textit{2}) \cdot p({\boldsymbol \theta}^\textit{2}\ \vert \ {\boldsymbol \alpha}^\textit{2}) \cdot p({\bf x}\ \vert \ {\bf z}^\textit{1}, {\bf z}^\textit{2}, {\boldsymbol \Phi}) \cdot p({\boldsymbol \Phi}\ \vert \ a_0, b_0) 
\end{align*}

第1項と第3項は
$$p({\bf z}\ \vert \ {\boldsymbol \theta}) = \prod_{n=1}^{N} \prod_{k=1}^{K} \theta_{k}^{\delta_{(z_n=k)}} = \prod_{k=1}^{K} \theta_{k}^{N_k}$$

第2項と第4項は
$$p({\boldsymbol \theta}\ \vert \ {\boldsymbol \alpha}) = \frac{\Gamma(\sum_{k=1}^{K} \alpha_k)}{\prod_{k=1}^{K} \Gamma(\alpha_k)} \prod_{k=1}^{K} \theta_k^{\alpha_k - 1} $$

第5項は
$$p({\bf x}\ \vert \ {\bf z}^\textit{1}, {\bf z}^\textit{2}, {\boldsymbol \Phi}) = \prod_{i=1}^{N^\textit{1}} \prod_{j=1}^{N^\textit{2}} \prod_{k=1}^{K} \prod_{\ell=1}^{L} \left[ \phi_{k\ell}^{x_{ij}} (1 - \phi_{k\ell})^{1 - x_{ij}} \right]^{\delta_{(z_i^\textit{1}=k, z_j^\textit{2}=\ell)}} = \prod_{k=1}^{K} \prod_{\ell=1}^{L} \phi_{k\ell}^{N_{k\ell}^{(+)}} (1 - \phi_{k\ell})^{N_{k\ell}^{(-)}}$$

第6項は
$$p({\boldsymbol \Phi}\ \vert \ a_0, b_0)= \prod_{k=1}^{K} \prod_{\ell=1}^{L} \frac{\Gamma(a_0+b_0)}{\Gamma(a_0)\Gamma(b_0)} \phi_{k\ell}^{a_0 - 1} (1 - \phi_{k\ell})^{b_0 - 1}$$

となります.
なお,$N_k = \sum_{n} \delta_{(z_n=k)}$,$N_{k\ell}^{(+)} = \sum_{i} \sum_{j} x_{ij} \delta_{(z_i^\textit{1}=k, z_j^\textit{2}=\ell)}$,$N_{k\ell}^{(-)} = \sum_{i} \sum_{j} (1 - x_{ij}) \delta_{(z_i^\textit{1}=k, z_j^\textit{2}=\ell)}$ と計算されます.

モデル推論

潜在クラスタ ${\bf z}$ の推論方法として collapsed Gibbs sampling(CGS)を紹介します.
collapsed とはモデルパラメータ $\{{\boldsymbol \theta}, {\boldsymbol \Phi}\}$ を積分消去(周辺化)することを意味します.

周辺同時確率は
$$p({\bf x}, {\bf z}^\textit{1}, {\bf z}^\textit{2}\ \vert \ {\boldsymbol \alpha}^\textit{1}, {\boldsymbol \alpha}^\textit{2}, a_0, b_0) = p({\bf z}^\textit{1}\ \vert \ {\boldsymbol \alpha}^\textit{1}) \cdot p({\bf z}^\textit{2}\ \vert \ {\boldsymbol \alpha}^\textit{2}) \cdot p({\bf x}\ \vert \ {\bf z}^\textit{1}, {\bf z}^\textit{2}, a_0, b_0)$$

であり,第1項と第2項は
$$p({\bf z}\ \vert \ {\boldsymbol \alpha}) = \int p({\bf z}\ \vert \ {\boldsymbol \theta}) \cdot p({\boldsymbol \theta}\ \vert \ {\boldsymbol \alpha}) d{\boldsymbol \theta} = \frac{\Gamma(\sum_{k=1}^{K} \alpha_k)}{\prod_{k=1}^{K} \Gamma(\alpha_k)} \frac{\prod_{k=1}^{K} \Gamma(N_k + \alpha_k)}{\Gamma(N + \sum_{k=1}^{K} \alpha_k)}$$

第3項は
$$p({\bf x}\ \vert \ {\bf z}^\textit{1}, {\bf z}^\textit{2}, a_0, b_0) = \int p({\bf x}\ \vert \ {\bf Z}, {\boldsymbol \Phi}) \cdot p({\boldsymbol \Phi}\ \vert \ a_0, b_0) d{\boldsymbol \Phi} = \prod_{k=1}^{K} \prod_{\ell=1}^{L} \frac{\Gamma(a_0+b_0)}{\Gamma(a_0)\Gamma(b_0)} \frac{\Gamma(a_0 + N_{k\ell}^{(+)}) \Gamma(b_0 + N_{k\ell}^{(-)})}{\Gamma(a_0 + b_0 + N_{k\ell})}$$

となります.

CGSでは第1ドメイン(行)のクラスタ $\{z_i^\textit{1}\}$ と第2ドメイン(列)のクラスタ $\{z_j^\textit{2}\}$ を繰り返しサンプリングすることで事後分布を推論します.

$\{z_i^\textit{1}\}$ のサンプリング式は次のようになります.

$$p(z_i^\textit{1}=k\ \vert \ {\bf x}, {\bf z}_{\backslash i}^\textit{1}, {\bf z}^\textit{2}, {\boldsymbol \alpha}^\textit{1}, {\boldsymbol \alpha}^\textit{2}, a _0, b _0) \propto p(z _i^\textit{1}=k\ \vert \ {\bf z} _{\backslash i}^\textit{1}, {\boldsymbol \alpha}^\textit{1}) \cdot p({\bf x} _{i:}\ \vert \ {\bf x} _{\backslash i:}, z _i^\textit{1}=k, {\bf z} _{\backslash i}^\textit{1}, {\bf z}^\textit{2}, a _0, b _0)$$

第1項は
$$p(z _i^\textit{1}=k\ \vert \ {\bf z} _{\backslash i}^\textit{1}, {\boldsymbol \alpha}^\textit{1}) = \frac{N _{k\backslash i}^\textit{1} + \alpha _k^\textit{1}}{N^\textit{1} - 1 + \sum _{k=1}^{K} \alpha _k^\textit{1}}$$

第2項は
$$p({\bf x} _{i:}\ \vert \ {\bf x} _{\backslash i:}, z _i^\textit{1}=k, {\bf z} _{\backslash i}^\textit{1}, {\bf z}^\textit{2}, a _0, b _0) = \prod _{\ell=1}^{L} \frac{\Gamma(a _0 + b _0 + N _{k\ell \backslash i:})}{\Gamma(a _0 + N _{k\ell \backslash i:}^{(+)}) \Gamma(b _0 + N _{k\ell \backslash i:}^{(-)})} \frac{\Gamma(a _0 + N _{k\ell \backslash i:}^{(+)} + N _{i\ell}^{(+)}) \Gamma(b _0 + N _{k\ell \backslash i:}^{(-)} + N _{i\ell}^{(-)})}{\Gamma(a _0 + b _0 + N _{k\ell \backslash i:} + N _\ell^\textit{2})}$$

なお,$N_{k\backslash n} = \sum_{n'\neq n} \delta_{(z_{n'}=k)}$,$N_{k\ell \backslash i:} = \sum_{i'\neq i} \sum_{j} \delta_{(z_{i'}^\textit{1}=k, z_j^\textit{2}=\ell)}$,$N_{i\ell}^{(+)} = \sum_{j} x_{ij} \delta_{(z_j^\textit{2}=\ell)}$,$N_{i\ell}^{(-)} = \sum_{j} (1 - x_{ij}) \delta_{(z_j^\textit{2}=\ell)}$ と計算されます.

$\{z_j^\textit{2}\}$ は上と対称なサンプリング式となります.

$$p(z _j^\textit{2}=\ell\ \vert \ {\bf x}, {\bf z}^\textit{1}, {\bf z} _{\backslash j}^\textit{2}, {\boldsymbol \alpha}^\textit{1}, {\boldsymbol \alpha}^\textit{2}, a _0, b _0) \propto p(z _j^\textit{2}=\ell\ \vert \ {\bf z} _{\backslash j}^\textit{2}, {\boldsymbol \alpha}^\textit{2}) \cdot p({\bf x} _{:j} \ \vert \ {\bf x} _{\backslash :j}, {\bf z}^\textit{1}, z _j^\textit{2}=\ell, {\bf z} _{\backslash j}^\textit{2}, a _0, b _0)$$

モデル実装

SBMのスクラッチ実装

以下は [関係データ行列x_, CGSのイテレーション数, クラスタ数($K$, $L$)] を引数としたときのSBMのスクラッチ実装です.

using StatsBase

const α1 = 0.5
const α2 = 0.5
const a0 = 0.5
const b0 = 0.5

# estimate model parameters with collapsed Gibbs sampling
function cgs(x_::Array, max_iter::Int, n_components::Tuple)
    N = size(x_)
    z_i  = zeros(Int, N[1])
    z_j  = zeros(Int, N[2])
    N_k  = zeros(Int, n_components[1])
    N_ℓ  = zeros(Int, n_components[2])
    N_kℓ = zeros(Int, n_components[1], n_components[2], 2)
    N_iℓ = zeros(Int, N[1], n_components[2], 2)
    N_jk = zeros(Int, N[2], n_components[1], 2)

    srand(0)
    @inbounds for iter in 1:max_iter

        # sampling z_i
        @inbounds for i in 1:N[1]
            x_i = x_[i,:]

            # remove x_i’s statistics
            if z_i[i] != 0
                k = z_i[i]
                N_k[k] -= 1
                @inbounds for (j, x) in enumerate(x_i)
                    N_jk[j, k, :] -= [1 - x, x]

                    if z_j[j] != 0
                         = z_j[j]
                        N_kℓ[k, , :] -= [1 - x, x]
                    end
                end
            end

            log_p_k = zeros(n_components[1])
            @inbounds for k in 1:n_components[1]
                log_p_k[k] = log(N_k[k] + α1)
                @inbounds for  in 1:n_components[2]
                    log_p_k[k] += lgamma(a0 + b0 + sum(N_kℓ[k, , :]))
                    log_p_k[k] -= lgamma(a0 + N_kℓ[k, , 2])
                    log_p_k[k] -= lgamma(b0 + N_kℓ[k, , 1])
                    log_p_k[k] += lgamma(a0 + N_kℓ[k, , 2] + N_iℓ[i, , 2])
                    log_p_k[k] += lgamma(b0 + N_kℓ[k, , 1] + N_iℓ[i, , 1])
                    log_p_k[k] -= lgamma(a0 + b0 + sum(N_kℓ[k, , :]) + N_ℓ[])
                end
            end

            # sample z_i after normalizing
            p_k = exp.(log_p_k - maximum(log_p_k))
            z_i[i] = sample(1:n_components[1], WeightVec(p_k / sum(p_k)))

            # add x_i’s statistics
            k = z_i[i]
            N_k[k] += 1
            @inbounds for (j, x) in enumerate(x_i)
                N_jk[j, k, :] += [1 - x, x]

                if z_j[j] != 0
                     = z_j[j]
                    N_kℓ[k, , :] += [1 - x, x]
                end
            end
        end

        # sampling z_j
        @inbounds for j in 1:N[2]
            x_j = x_[:,j]

            # remove transpose(x)_j’s statistics
            if z_j[j] != 0
                 = z_j[j]
                N_ℓ[] -= 1
                @inbounds for (i, x) in enumerate(x_j)
                    N_iℓ[i, , :] -= [1 - x, x]

                    if z_i[i] != 0
                        k = z_i[i]
                        N_kℓ[k, , :] -= [1 - x, x]
                    end
                end
            end

            log_p_ℓ = zeros(n_components[2])
            @inbounds for  in 1:n_components[2]
                log_p_ℓ[] = log(N_ℓ[] + α2)
                @inbounds for k in 1:n_components[1]
                    log_p_ℓ[] += lgamma(a0 + b0 + sum(N_kℓ[k, , :]))
                    log_p_ℓ[] -= lgamma(a0 + N_kℓ[k, , 2])
                    log_p_ℓ[] -= lgamma(b0 + N_kℓ[k, , 1])
                    log_p_ℓ[] += lgamma(a0 + N_kℓ[k, , 2] + N_jk[j, k, 2])
                    log_p_ℓ[] += lgamma(b0 + N_kℓ[k, , 1] + N_jk[j, k, 1])
                    log_p_ℓ[] -= lgamma(a0 + b0 + sum(N_kℓ[k, , :]) + N_k[k])
                end
            end

            # sample z_j after normalizing
            p_ℓ = exp.(log_p_ℓ - maximum(log_p_ℓ))
            z_j[j] = sample(1:n_components[2], WeightVec(p_ℓ / sum(p_ℓ)))

            # add transpose(x)_j’s statistics
             = z_j[j]
            N_ℓ[] += 1
            @inbounds for (i, x) in enumerate(x_j)
                N_iℓ[i, , :] += [1 - x, x]

                if z_i[i] != 0
                    k = z_i[i]
                    N_kℓ[k, , :] += [1 - x, x]
                end
            end
        end
    end
end

上図の関係データ行列でモデルを実行すると

julia> x_
15×10 Array{Int64,2}:
 1  0  0  1  1  1  0  0  0  0
 0  1  1  0  1  0  1  0  0  1
 1  1  1  0  1  1  0  0  0  1
 1  0  0  1  0  0  1  1  1  0
 1  0  0  0  1  1  0  0  0  0
 1  0  0  0  1  1  0  0  0  1
 0  0  1  1  0  0  0  1  1  0
 0  0  0  0  0  0  0  1  1  0
 1  1  0  0  1  0  1  0  0  0
 0  1  0  0  1  1  0  0  0  0
 1  1  1  0  1  1  1  0  0  1
 1  0  0  0  1  1  0  0  1  0
 0  0  0  1  0  0  0  1  1  0
 1  0  0  0  1  0  0  0  0  0
 0  0  0  1  0  0  0  1  1  0

julia> fit(x_, max_iter=100, n_components=(3,3))

$\{z_i^\textit{1}\} = [3, 1, 1, 2, 3, 3, 2, 2, 1, 3, 1, 3, 2, 3, 2]$
$\{z_j^\textit{2}\} = [2, 1, 1, 3, 2, 2, 1, 3, 3, 1]$

となり,ソートすると右図のブロック構造にクラスタリングできています.

IRMのスクラッチ実装

SBMではクラスタ数 $(K,L)$ を事前に決める必要がありますが,SBMのノンパラベイズ拡張であるIRM(infinite relational model; 無限関係モデル)では適切なクラスタ数を推定することができます.

IRMの生成過程は以下のようになります.

\begin{align}
{\boldsymbol \theta}^\textit{1} &\sim {\rm SBP}(\alpha^\textit{1}) \\
{\boldsymbol \theta}^\textit{2} &\sim {\rm SBP}(\alpha^\textit{2}) \\
\{\{ \phi_{k\ell} \}_{\ell=1}^{\infty} \}_{k=1}^{\infty} &\sim {\rm Beta}(a_0, b_0) \\
\{ z_i^\textit{1} \}_{i=1}^{N^\textit{1}} &\sim {\rm Categorical}({\boldsymbol \theta}^\textit{1}) \\
\{ z_j^\textit{2} \}_{j=1}^{N^\textit{2}} &\sim {\rm Categorical}({\boldsymbol \theta}^\textit{2}) \\
\{\{ x_{ij} \}_{j=1}^{N^\textit{2}} \}_{i=1}^{N^\textit{1}} &\sim {\rm Bernoulli}(\phi_{z_i^\textit{1} z_j^\textit{2}})
\end{align}

SBP(stick-breaking process; 棒折り過程)は無限次元の離散分布から生成されることを意味します.
以下はIRMのスクラッチ実装です.

function cgs(x_::Array, max_iter::Int)
    N = size(x_)
    z_i  = zeros(Int, N[1])
    z_j  = zeros(Int, N[2])
    N_k  = zeros(Int, 0)
    N_ℓ  = zeros(Int, 0)
    N_kℓ = zeros(Int, 0, 0, 2)
    N_iℓ = zeros(Int, N[1], 0, 2)
    N_jk = zeros(Int, N[2], 0, 2)
    n_components = [0, 0]

    srand(0)
    @inbounds for iter in 1:max_iter

        # sampling z_i
        @inbounds for i in 1:N[1]
            x_i = x_[i,:]

            # remove x_i’s statistics
            if z_i[i] != 0
                k = z_i[i]
                N_k[k] -= 1
                @inbounds for (j, x) in enumerate(x_i)
                    N_jk[j, k, :] -= [1 - x, x]

                    if z_j[j] != 0
                         = z_j[j]
                        N_kℓ[k, , :] -= [1 - x, x]
                    end
                end

                # if any component is empty, remove it and decrease n_components
                if N_k[k] == 0
                    n_components[1] -= 1
                    deleteat!(N_k, k)
                    N_jk = N_jk[:, setdiff(1:end, k), :]
                    N_kℓ = N_kℓ[setdiff(1:end, k), :, :]
                    @inbounds for (tmp_i, tmp_z) in enumerate(z_i)
                        if tmp_z > k
                            z_i[tmp_i] -= 1
                        end
                    end
                end
            end

            K       = n_components[1] + 1
            log_p_k = zeros(K)

            # existing components
            @inbounds for k in 1:n_components[1]
                log_p_k[k] = log(N_k[k])
                @inbounds for  in 1:n_components[2]
                    log_p_k[k] += lgamma(a0 + b0 + sum(N_kℓ[k, , :]))
                    log_p_k[k] -= lgamma(a0 + N_kℓ[k, , 2])
                    log_p_k[k] -= lgamma(b0 + N_kℓ[k, , 1])
                    log_p_k[k] += lgamma(a0 + N_kℓ[k, , 2] + N_iℓ[i, , 2])
                    log_p_k[k] += lgamma(b0 + N_kℓ[k, , 1] + N_iℓ[i, , 1])
                    log_p_k[k] -= lgamma(a0 + b0 + sum(N_kℓ[k, , :]) + N_ℓ[])
                end
            end

            # a new component
            log_p_k[K]  = log(α1)
            @inbounds for  in 1:n_components[2]
                log_p_k[K] += lgamma(a0 + b0)
                log_p_k[K] -= lgamma(a0)
                log_p_k[K] -= lgamma(b0)
                log_p_k[K] += lgamma(a0 + N_iℓ[i, , 2])
                log_p_k[K] += lgamma(b0 + N_iℓ[i, , 1])
                log_p_k[K] -= lgamma(a0 + b0 + N_ℓ[])
            end

            # sample z_i after normalizing
            p_k = exp.(log_p_k - maximum(log_p_k))
            z_i[i] = sample(1:K, WeightVec(p_k / sum(p_k)))

            if z_i[i] == K
                n_components[1] += 1
                push!(N_k, 0)
                N_jk = cat(2, N_jk, zeros(Int, N[2], 1, 2))
                N_kℓ = cat(1, N_kℓ, zeros(Int, 1, n_components[2], 2))
            end

            # add x_i’s statistics
            k = z_i[i]
            N_k[k] += 1
            @inbounds for (j, x) in enumerate(x_i)
                N_jk[j, k, :] += [1 - x, x]

                if z_j[j] != 0
                     = z_j[j]
                    N_kℓ[k, , :] += [1 - x, x]
                end
            end
        end

        # sampling z_j
        @inbounds for j in 1:N[2]
            x_j = x_[:,j]

            # remove transpose(x)_j’s statistics
            if z_j[j] != 0
                 = z_j[j]
                N_ℓ[] -= 1
                @inbounds for (i, x) in enumerate(x_j)
                    N_iℓ[i, , :] -= [1 - x, x]

                    if z_i[i] != 0
                        k = z_i[i]
                        N_kℓ[k, , :] -= [1 - x, x]
                    end
                end

                # if any component is empty, remove it and decrease n_components
                if N_ℓ[] == 0
                    n_components[2] -= 1
                    deleteat!(N_ℓ, )
                    N_iℓ = N_iℓ[:, setdiff(1:end, ), :]
                    N_kℓ = N_kℓ[:, setdiff(1:end, ), :]
                    @inbounds for (tmp_j, tmp_z) in enumerate(z_j)
                        if tmp_z > 
                            z_j[tmp_j] -= 1
                        end
                    end
                end
            end

            L       = n_components[2] + 1
            log_p_ℓ = zeros(L)

            # existing components
            @inbounds for  in 1:n_components[2]
                log_p_ℓ[] = log(N_ℓ[] + α2)
                @inbounds for k in 1:n_components[1]
                    log_p_ℓ[] += lgamma(a0 + b0 + sum(N_kℓ[k, , :]))
                    log_p_ℓ[] -= lgamma(a0 + N_kℓ[k, , 2])
                    log_p_ℓ[] -= lgamma(b0 + N_kℓ[k, , 1])
                    log_p_ℓ[] += lgamma(a0 + N_kℓ[k, , 2] + N_jk[j, k, 2])
                    log_p_ℓ[] += lgamma(b0 + N_kℓ[k, , 1] + N_jk[j, k, 1])
                    log_p_ℓ[] -= lgamma(a0 + b0 + sum(N_kℓ[k, , :]) + N_k[k])
                end
            end

            # a new component
            log_p_ℓ[L]  = log(α2)
            @inbounds for k in 1:n_components[1]
                log_p_ℓ[L] += lgamma(a0 + b0)
                log_p_ℓ[L] -= lgamma(a0)
                log_p_ℓ[L] -= lgamma(b0)
                log_p_ℓ[L] += lgamma(a0 + N_jk[j, k, 2])
                log_p_ℓ[L] += lgamma(b0 + N_jk[j, k, 1])
                log_p_ℓ[L] -= lgamma(a0 + b0 + N_k[k])
            end

            # sample z_j after normalizing
            p_ℓ = exp.(log_p_ℓ - maximum(log_p_ℓ))
            z_j[j] = sample(1:L, WeightVec(p_ℓ / sum(p_ℓ)))

            if z_j[j] == L
                n_components[2] += 1
                push!(N_ℓ, 0)
                N_iℓ = cat(2, N_iℓ, zeros(Int, N[1], 1, 2))
                N_kℓ = cat(2, N_kℓ, zeros(Int, n_components[1], 1, 2))
            end

            # add transpose(x)_j’s statistics
             = z_j[j]
            N_ℓ[] += 1
            @inbounds for (i, x) in enumerate(x_j)
                N_iℓ[i, , :] += [1 - x, x]

                if z_i[i] != 0
                    k = z_i[i]
                    N_kℓ[k, , :] += [1 - x, x]
                end
            end
        end
    end
end

こちらも上と同様の関係データ行列でモデルを実行すると

$\{z_i^\textit{1}\} = [3, 1, 1, 2, 3, 3, 2, 2, 1, 3, 1, 3, 2, 3, 2]$
$\{z_j^\textit{2}\} = [3, 2, 2, 1, 3, 3, 2, 1, 1, 2]$

となり,$(K,L)=(3,3)$ のブロック構造にクラスタリングすることができました.

実用的には事前にクラスタ数がわかっている場面は多くありません.
そんなときにはクラスタ数を自動決定できるノンパラベイズモデルが非常に有用な手段となるでしょう.

まとめ

確率的ブロックモデルについて紹介しました.
スクラッチ実装したSBMとIRMは https://github.com/yusaku-i/Bayesian_mixtures.jl/tree/master/src/SBM で公開しています.

最後にJulia×データ分析に関する執筆依頼お待ちしております.

Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.