Julia Advent Calendar 2017の16日目の記事です.
内容は「GMM (Gaussian mixture model; 混合ガウスモデル) を実装してパッケージ公開した話」です.

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

Juliaはまだまだ日本語リファレンスが乏しいため,Pythonに比べるとググり力が必要になります.
そこで今回は,Julia初学者のためにモデルの実装からパッケージとしてGitHubに公開するまでの一連の流れについて紹介したいと思います.

環境設定

まずはJuliaの環境設定を行います.
インストールがまだの方はここからOS環境に応じて必要なファイルをダウンロードしてください.
macをお使いの方は以下のbrewコマンドでも最新verをインストールできます.

$ brew cask install julia

インストールが完了したらJuliaを起動しましょう.

スクリーンショット 2017-12-06 15.49.25.png

次にパッケージ作成のためにPkgDev.jlを追加します.

julia> Pkg.add("PkgDev")

追加が終わったら以下のコマンドを入力します.

julia> using PkgDev

julia> PkgDev.generate("hoge", "MIT")

すると

スクリーンショット 2017-12-06 16.04.49.png

のようにパッケージ公開に必要なファイルが生成されます.
生成されたフォルダの中身は

スクリーンショット 2017-12-06 16.12.27.png

のような構成になっています.

モデル実装

モデルの実装はhoge/src/hoge.jlで行います.
まずは"Hello World!"を出力するように実装します.

スクリーンショット 2017-12-06 16.42.48.png

以下のコマンドでhogeモジュールを読み出して"Hello World!"を出力します.

スクリーンショット 2017-12-06 16.44.06.png

無事に"Hello World!"を出力できました.
あとはモデルの実装をひたすら進めるだけです.
以下はGMMのスクラッチ実装です.

using StatsBase

const α  = 1.0
const m0 = 0.0
const ξ0 = 1.0
const a0 = 1.0
const b0 = 1.0

const μ0 = m0
const λ0 = (ξ0 * a0) / ((1 + ξ0) * b0)
const ν0 = 2 * a0

# estimate model parameters with collapsed Gibbs sampling
function fit(x_::Vector; max_iter::Int=100, n_components=nothing)
    @assert max_iter > 0

    # finite Gaussian mixture model
    if n_components != nothing
        @assert n_components > 0

        finite_cgs(x_, max_iter, n_components)

    # infinite Gaussian mixture model
    else
        infinite_cgs(x_, max_iter)
    end
end

# compute the per-sample average log-likelihood of the given data
function score(x_::Vector)
    @assert mode  ["fin", "inf"]

    ret = 0.0
    @inbounds for x in x_
        ret += log(dot(θ_, [normal_pdf(x, μ_[k], λ_[k]) for k in 1:n_components]))
    end
    ret /= length(x_)

    return ret
end

# predict the labels for the given data using trained model
function predict(x_::Vector)
    @assert mode  ["fin", "inf"]

    ret = Vector{Int}()
    local denom = sum(D_k) + α
    @inbounds for x in x_
        local log_t_ = [t_logpdf(x, μ_[k], λ_[k], ν_[k]) for k in 1:n_components]

        if mode == "fin"
            push!(ret, indmax(θ_ .* exp.(log_t_)))
        elseif mode == "inf"
            let
                local θ_ = [D_k[k] / denom for k in 1:n_components]
                push!(θ_, α / denom)
                push!(log_t_, t_logpdf(x, μ0, λ0, ν0))
                push!(ret, indmax(θ_ .* exp.(log_t_)))
            end
        end
    end

    return ret
end

function finite_cgs(x_::Vector, max_iter::Int, n_components::Int)
    local z_d = zeros(Int, length(x_))
    global D_k = zeros(Int, n_components)
    local C_k = zeros(n_components)
    local S_k = zeros(n_components)

    srand(0)
    @inbounds for iter in 1:max_iter
        @inbounds for (d, x) in enumerate(x_)

            # remove x_d’s statistics
            if z_d[d] != 0
                D_k[z_d[d]] -= 1
                C_k[z_d[d]] -= x
                S_k[z_d[d]] -= x^2
            end

            local log_p_z = zeros(n_components)
            @inbounds for k in 1:n_components
                log_p_z[k] = log(D_k[k] + α)

                ξ = D_k[k] + ξ0
                m = (C_k[k] + ξ0 * m0) / ξ
                a = D_k[k] / 2 + a0
                b = (S_k[k] + ξ0 * m0^2 - ξ * m^2) / 2 + b0

                μ = m
                λ = (ξ * a) / ((1 + ξ) * b)
                ν = 2 * a

                log_p_z[k] += t_logpdf(x, μ, λ, ν)
            end

            # sample z_d after normalizing
            local p_z = exp.(log_p_z - maximum(log_p_z))
            z_d[d] = sample(1:n_components, WeightVec(p_z / sum(p_z)))

            # add x_d’s statistics
            D_k[z_d[d]] += 1
            C_k[z_d[d]] += x
            S_k[z_d[d]] += x^2
        end
    end

    posteriori_estimation(n_components, D_k, C_k, S_k, "Dir")
end

function infinite_cgs(x_::Vector, max_iter::Int)
    local z_d = zeros(Int, length(x_))
    global D_k = zeros(Int, 0)
    local C_k = zeros(0)
    local S_k = zeros(0)
    local n_components = 0

    srand(0)
    @inbounds for iter in 1:max_iter
        @inbounds for (d, x) in enumerate(x_)

            # remove x_d’s statistics
            if z_d[d] != 0
                D_k[z_d[d]] -= 1
                C_k[z_d[d]] -= x
                S_k[z_d[d]] -= x^2

                # if any component is empty, remove it and decrease n_components
                if D_k[z_d[d]] == 0
                    n_components -= 1
                    deleteat!(D_k, z_d[d])
                    deleteat!(C_k, z_d[d])
                    deleteat!(S_k, z_d[d])
                    local z = z_d[d]
                    @inbounds for (tmp_d, tmp_z) in enumerate(z_d)
                        if tmp_z > z
                            z_d[tmp_d] -= 1
                        end
                    end
                end
            end

            local K = n_components + 1
            local log_p_z = zeros(K)

            # existing components
            @inbounds for k in 1:n_components
                log_p_z[k] = log(D_k[k] + α)

                ξ = D_k[k] + ξ0
                m = (C_k[k] + ξ0 * m0) / ξ
                a = D_k[k] / 2 + a0
                b = (S_k[k] + ξ0 * m0^2 - ξ * m^2) / 2 + b0

                μ = m
                λ = (ξ * a) / ((1 + ξ) * b)
                ν = 2 * a

                log_p_z[k] += t_logpdf(x, μ, λ, ν)
            end

            # a new component
            log_p_z[K] = log(α)
            log_p_z[K] += t_logpdf(x, μ0, λ0, ν0)

            # sample z_d after normalizing
            local p_z = exp.(log_p_z - maximum(log_p_z))
            z_d[d] = sample(1:K, WeightVec(p_z / sum(p_z)))

            if z_d[d] == K
                n_components += 1
                push!(D_k, 0)
                push!(C_k, 0.0)
                push!(S_k, 0.0)
            end

            # add x_d’s statistics
            D_k[z_d[d]] += 1
            C_k[z_d[d]] += x
            S_k[z_d[d]] += x^2
        end
    end

    posteriori_estimation(n_components, D_k, C_k, S_k, "SBP")
end

function posteriori_estimation(K::Int, D_k::Vector, C_k::Vector, S_k::Vector, prior::String)
    @assert prior  ["Dir", "SBP"]

    global θ_ = zeros(K)
    global μ_ = zeros(K)
    global λ_ = zeros(K)
    global ν_ = zeros(K)
    global n_components = K

    global mode = ""
    if prior == "Dir"
        mode = "fin"
    elseif prior == "SBP"
        mode = "inf"
    end

    local D = sum(D_k)
    @inbounds for k in 1:K
        if prior == "Dir"
            θ_[k] = (D_k[k] + α) / (D + α * K)
        elseif prior == "SBP"
            θ_[k] = D_k[k] / D
        end

        ξ = D_k[k] + ξ0
        m = (C_k[k] + ξ0 * m0) / ξ
        a = D_k[k] / 2 + a0
        b = (S_k[k] + ξ0 * m0^2 - ξ * m^2) / 2 + b0

        μ_[k] = m
        λ_[k] = (ξ * a) / ((1 + ξ) * b)
        ν_[k] = 2 * a
    end
end

function t_logpdf(x::Real, μ::Float64, λ::Float64, ν::Float64)
    # ret = gamma((ν + 1) / 2) / gamma(ν / 2) * (λ / (π * ν))^(0.5) * (1 + λ / ν * (x - μ)^2)^(- (ν + 1) / 2)
    ret = lgamma((ν + 1) / 2) - lgamma(ν / 2) + 0.5 * log(λ / (π * ν)) - (ν + 1) / 2 * log(1 + λ / ν * (x - μ)^2)
    return ret
end

function normal_pdf(x::Real, μ::Float64, λ::Float64)
    ret = (λ / (2 * π))^(0.5) * exp(- λ / 2 * (x - μ)^2)
    return ret
end

なお,今回は1次元のデータ(スカラー値)に対するGMMです.

GMMは一言で言えば「ガウス分布の線形重ね合わせで表されるモデル」になります.
構成要素 $K$ 個のGMMは以下の式になります.

p(x) = \sum_{k=1}^{K} \theta_k \mathcal{N} (x\ \vert \ \mu_k, \lambda_k^{-1}) 

ここで,$\theta_k$ はクラス混合比を意味するパラメータ, $\mu_k$ と $\lambda_k$ はそれぞれガウス分布の平均と精度 (分散の逆数) を意味するパラメータになります.

GMMを詳しく知りたい方は
- パターン認識と機械学習
- ベイズ推論による機械学習入門
などを参考にしてください.

パッケージ公開

ここまで来たらあとは実装したモデルをパッケージとしてGitHubに公開するだけです.

まずはhoge.jlという名前でレポジトリを作成してください.
GitHubアカウントは以下のコマンドで確認できます.

julia> PkgDev.config()

レポジトリ作成が完了したら以下のコマンドで実装したモデルをpushします.

git add src/hoge.jl
git commit -m "first commit"
git remote add origin git@github.com:$YOUR_GITHUB_ACCOUNT/hoge.jl.git
git push -u origin master

すると

スクリーンショット 2017-12-06 19.01.02.png

これで世界中のプログラマーがhoge.jlをインストールできるようになりました.
ローカル環境へのインストール方法は以下のコマンドになります.

julia> Pkg.clone("git@github.com:$YOUR_GITHUB_ACCOUNT/hoge.jl.git")

ちなみにJuliaの公式レポジトリに自家製パッケージを公開したい場合は METADATA.jlプルリクを送る いろいろゴニョゴニョする 必要があります.
無事にマージされると

julia> Pkg.add("hoge")

でインストール可能になります.

まとめ

「GMMを実装してパッケージ公開した話」を紹介しました.
と言ってもほとんどJuliaのパッケージ公開のための手順書のようになってしまいました.

ちなみにスクラッチ実装したGMMは
https://github.com/yusaku-i/UnivariateGaussianMixtureModel.jl
(指摘を受けて修正)https://github.com/yusaku-i/Bayesian_mixtures.jl/tree/master/src/uGMM で公開しています
また、多次元GMMも https://github.com/yusaku-i/Bayesian_mixtures.jl/tree/master/src/mGMM で公開しています、ぜひ触ってみてください.

まだまだマイナーな言語の印象のあるJuliaですが,今後データ分析界隈では最も使われる言語になること間違いなしです!
最後にJulia×データ分析に関する執筆依頼お待ちしております笑

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