[Package] Optim.jlを使って最尤推定 | コイン・多項ロジット
まえがき
最尤推定法 (maximum likelihood method) は機械学習における基本的なアプローチとして頻出します.細かい話は置いておくとして,観測したデータから尤度関数
最尤推定 | コインの例
例題
よくある例題です.コインを振って表が出る確率を
となります.最尤推定法ではこれを最大化する
Julia Plots.jlを用いた確認
ここでは
using Plots
gr()
Lterm(x; N=6, r=5) = x ^ r * (1 - x) ^ (N - r);
N, r = 6, 5
xseq = 0:0.01:1.0;
yseq = Lterm.(xseq);
f = plot(size=(400, 300), legend=:outertopright)
plot!(f, xseq, yseq, color=:tomato, lw=3, label="likelihood")
plot!(f, [r/N, r/N], [0, Lterm(r/N)], label="est")
f
Julia Optim.jlを利用した推定
今回は閉じた式
以下が実装です.対数尤度の計算で
using Optim
using Plots
gr()
# 対数尤度関数 Optim.jl 用
LL(q; N=6, r=5) = -(log(N) + r * log(q[1]) + (N - r) *log(1 - q[1]))
# Optim.jlで最適値を求める
res = optimize(LL, [0.0001], [1.0], [0.5])
xm = Optim.minimizer(res)[1];
# 描画
xseq = 0:0.01:1.0;
yseq = LL.(xseq);
f = plot(size=(400, 300), legend=:outertopright)
plot!(f, xseq, yseq, color="tomato", lw=3, label="log-likelihood")
plot!(f, [xm, xm], [0.0, LL(xm)], lw=5, label="min")
f
最尤推定 | 意思決定における多項ロジットの例
離散選択モデル (いくつかの候補が離散的に与えられ,そのうち1つを選ぶ) は意思決定の基本的なモデルとして経済学 (計量経済学) などを中心に広く使われています. 理論的な背景は書籍やインターネット上の文書にお願いするとし,ここでは簡単な例だけを考えます.数式の詳細は置いておくとして,
で計算するというもので,ニューラルネットワークなどを日々扱っている人には馴染みの深い計算かと思います.この場合の観測データとしては,
例題
ここで例題を用いて具体的な例を使って計算をしてみます.例題は参考文献から取ってきています.ユーザが適当にいたとして選択肢として
の式で書かれていると仮定します.このときいろんな人の選好が混じっており,同じ
参考文献から拝借してきたデータは次の表です.
計算
上のモデルにおいて
今データに値
になるとのことです.
Julia Optim.jlを利用した推定
雑に実装します (真似しないでください).今回は
using Optim
using Plots
gr()
# データ
data = [
20 5 210 50 10 120;
20 10 210 100 30 90;
30 15 210 100 10 60;
30 20 210 125 10 60;
45 30 210 150 30 90;
45 45 210 250 60 70;
60 60 420 420 20 20;
60 30 420 1100 25 10;
60 45 420 1100 30 5
]
# 効用と確率
V1n(x, α, β) = α * x[1] + β * x[2];
V2n(x, α, β, γ) = α * x[1] + β * x[2] + γ;
P1n(x, α, β, γ) = exp(α * x[1] + β * x[2]) / (exp(α * x[1] + β * x[2]) + exp(α * x[3] + β * x[4] + γ));
# 対数尤度 (最小化問題にするために-符号)
function LL(param)
α, β, γ = param
ll = 0.0
for row in eachrow(data)
p1n = P1n([row[1] row[3] row[2] row[4]], α, β, γ)
p2n = 1 - p1n
l1 = row[5] * log(p1n)
l2 = row[6] * log(p2n)
ll += l1 + l2
end
-ll
end
# Optim.jlで最小化
res = optimize(eval, [0.0, 0.0, 0.0], autodiff=:forward, LBFGS())
# 出力
α, β, γ = Optim.minimizer(res)
println(α)
println(β)
println(γ)
println(LL(Optim.minimizer(res)))
特に終了条件などをいじらず適当に使ったOptim.jlの計算によれば
となりました.Excelの結果と少し異なりますが,だいたいうまく動いたと思って良さそうですね.
参考文献
- Optim.jlについては公式のドキュメントを参照しました: https://julianlsolvers.github.io/Optim.jl/stable/#
- 多項ロジットの数値例についてはこちらのスライドから拝借しました: http://bin.t.u-tokyo.ac.jp/startup14/file/2-1.pdf
- 他には世の中の説明などを参考しました: https://h-memo.com/logit-model/
Discussion