JuliaでEMアルゴリズムの関数を作る

なぜやったか

• 仕事上混合正規分布のパラメータ推定への理解が必要となってしまいました。
  • パラメータ推定によく用いられているEMアルゴリズムを理解したい
  • せっかくなので正規分布だけでなくポアソン分布などのケースも考えました。

参考文献

1.は、EMアルゴリズムの理解のために参考としたものです。シミュレーションの定数も、ある程度参考としています。ありがとうございます。

1. https://qiita.com/gorogoro_gorori/items/fdecda8e237aef74e3d0

やったこと

1. 参考文献の1.を読みました。

2. Juliaで一次元の混合正規分布、混合ポアソン分布、混合指数分布の推定用の関数を作成しました。実際に、混合ポアソン分布の推定結果をこの記事で整理しています。

コード

# パッケージの呼び出し
using Distributions
using LinearAlgebra

# 混合正規分布のパラメータ推定を行う関数
function em_gaussian_1dim(x, iter_num, pis, mus, sigmas)
    N = length(x)
    for ite in 1:iter_num
        piNs = [p * pdf.(Normal(m, s), x) for (p, m, s) in zip(pis, mus, sigmas)];
        qns = [pn ./ sum(piNs) for pn in piNs];
        pis = [sum(qn) / N for qn in qns];
        mus = [dot(qn, x) / (N * p) for (p, qn) in zip(pis, qns)];
        sigmas = [sqrt(sum(qn .* (x .- m).^2) / (N * p)) for (m, p, qn) in zip(mus, pis, qns)];
    end
    return pis, mus, sigmas
end

# 混合指数分布用のパラメータ推定を行う関数
function em_exponential(x, iter_num, pis, lambdas)
    N = length(x)
    for ite in 1:iter_num
        piNs = [p * pdf.(Exponential(lambda), x) for (p, lambda) in zip(pis, lambdas)];
        qns = [pn ./ sum(piNs) for pn in piNs];
        pis = [sum(qn) / N for qn in qns];
        lambdas = [dot(qn, x) / (N * p) for (p, qn) in zip(pis, qns)];
    end
    return pis, lambdas
end

# 混合ポアソン分布のパラメータ推定を行う関数
function em_poisson(x, iter_num, pis, lambdas)
    N = length(x)
    for ite in 1:iter_num
        piNs = [p * pdf.(Poisson(lambda), x) for (p, lambda) in zip(pis, lambdas)];
        qns = [pn ./ sum(piNs) for pn in piNs];
        pis = [sum(qn) / N for qn in qns];
        lambdas = [dot(qn, x) / (N * p) for (p, qn) in zip(pis, qns)];
    end
    return pis, lambdas
end

# 乱数生成
y0 = rand(DiscreteUniform(1,3), 100000);
x = reshape([rand(Poisson(y^2), 1)[1] for y in y0], (100000, 1));
iter_num = 1000;
pis = [1/2 1/3 1/6];
lambdas = [1 5 10];

# 実際の推定
em_poisson(x, iter_num, pis, lambdas)
# 出力結果
# 実行時間は30秒ほど、およそパラメータは推定できているっぽいです。
# 推定結果をアニメーションにするなどの方法を知りたいですがまた今度。
# ([0.3369903366831118 0.33647632001950306 0.3265333432973852], 
# [0.9977504488001773 4.072583397747915 9.044037695711244])

感想

• ある程度推定結果は正しいと感じています。
• 正規分布を多次元に拡張した場合や、実際のデータを使ったモデルの推定も記事にしたいです。

次にやる、この記事に関連すること

• 実際のデータを使ったモデルの推定の記事化。