確率的プログラミングモデルのGenのtutorialを触ってみた
Genとは
初心者向け汎用AIプログラミングシステムとの記事があって気になったのですが、数学の知識がなくても使えるとこことで普通のプログラマーでも扱えることを意識しているのかと思いました。Genの論文ではどのようなモデルでも表現できる『汎用の』モデリング言語を提供しているシステムもあるが、そういったシステムは使いものにならないほど収束の遅い推論アルゴリズムしかサポートしていないとのことで、Genでは汎用目的でもその問題を改善しようとしているらしいです。
環境準備
公式サイトにJulia package managerを使った方法とDockerを使った方法での環境構築が載っています。今回はJulia package managerを使って確認したいと思います。
juliaインストール(mac)
Genをインストールする前にまずjuliaを使えるようにする必要があります。ここからダウンロードすることができますが、macの場合は以下のコマンドで簡単にインストールできます。
brew cask install julia
Jupyter Notebookのカーネルにjuliaを追加する
Jupyter Notebookで動作確認したいのでカーネルを追加します。julia起動後に]を入力するとプロンプトがjulia>から(vx.x) pkg>に変わるので、それからadd IJuliaを実行するだけでJupyter Notebook上でjuliaが実行できるようになります。
$ julia
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.1.1 (2019-05-16)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
(v1.1) pkg>
(v1.1) pkg> add IJulia
Genインストール
それからGenをインストールします。
(v1.1) pkg> add https://github.com/probcomp/Gen
tutorialを試してみる
jupyter-notebook上でtutorialを試してみます。
using Gen using PyPlot
自分の環境だとusing PyPlotで以下のエラーが出たので、メッセージに従いimport Pkg; Pkg.add("PyPlot")を実行しておきました。
ArgumentError: Package PyPlot not found in current path:
- Run `import Pkg; Pkg.add("PyPlot")` to install the PyPlot package.
Stacktrace:
[1] require(::Module, ::Symbol) at ./loading.jl:823
[2] top-level scope at In[2]:1
末尾にセミコロンをつけなくても正しく動きますが、確認したい結果かどうかを区別したいので、確認する必要のないものについては末尾にセミコロンをつけておいて確認したいものはprintln関数を呼ぶなどします。
確率関数モデルを書く
Genで確率モデルを書く場合は組み込みDSLを使いデータ生成関数を奇術します。関数に@genがついていますが、これはこの関数がモデル化されたデータ生成プロセスということを表しています。
以下の関数ではxy平面上の線形モデルの確率モデルを表していて、一連のxのデータのみを与えただけで平面内の線を無造作に選び与えられたxに対してのy座標が線の近くになるように選ばれるらしい。
@gen function line_model(xs::Vector{Float64})
n = length(xs)
# We begin by sampling a slope and intercept for the line.
# Before we have seen the data, we don't know the values of
# these parameters, so we treat them as random choices. The
# distributions they are drawn from represent our prior beliefs
# about the parameters: in this case, that neither the slope nor the
# intercept will be more than a couple points away from 0.
slope = @trace(normal(0, 1), :slope)
intercept = @trace(normal(0, 2), :intercept)
# Given the slope and intercept, we can sample y coordinates
# for each of the x coordinates in our input vector.
for (i, x) in enumerate(xs)
@trace(normal(slope * x + intercept, 0.1), (:y, i))
end
# The return value of the model is often not particularly important,
# Here, we simply return n, the number of points.
return n
end;
slope = @trace(normal(0, 1), :slope)とintercept = @trace(normal(0, 2), :intercept)はそれぞれ平均0で標準偏差1の正規分布と平均0で標準偏差2の正規分布を表していて、@trace(normal(slope * x + intercept, 0.1), (:y, i))の部分ではそれぞれを線形モデルの傾き、切片としたものが平均で標準偏差0.１とした正規分布に従いyのデータを生成しているようです。
関数の戻り値自体はlength(xs)になっているので、以下を実行すると配列のサイズである11が表示される。
xs = [-5., -4., -3., -.2, -1., 0., 1., 2., 3., 4., 5.]; n = line_model(xs) println(n) 11
slope = @trace(normal(0, 1), :slope)ではtrace関数内の:slopeはデータ生成関数のトレース対象になるようですが生成した値を直接参照したい場合は変数に代入しておく必要があるようです。直後にprintln(slope)を実行させると関数呼び出しごとに値が変わるのが確認できます。trace関数の第二期引数で渡している値はGen.simulateを実行することで確認ができます。
trace = Gen.simulate(line_model, (xs,)); println(trace)
生成した値のみを確認したいのであればGen.get_choicesを使います。
println(Gen.get_choices(trace))
生成した値から特定のものを取得したい場合は以下のようになります。
choices = Gen.get_choices(trace) println(choices[:slope])
データ生成関数に渡した引数はGen.get_argsで確認できます。
Gen.get_args(trace)
データ生成関数の戻り値はGen.get_retvalで取得できます。
println(Gen.get_retval(trace));
データの確認はPython同様Matplotlibなど使います。今回は線形モデルを先に決めた上でyの座標を決めていたので確実に線の近くに点がつくのですが、
function render_trace(trace; show_data=true)
# Pull out xs from the trace
xs = get_args(trace)[1]
xmin = minimum(xs)
xmax = maximum(xs)
if show_data
ys = [trace[(:y, i)] for i=1:length(xs)]
# Plot the data set
scatter(xs, ys, c="black")
end
# Pull out slope and intercept from the trace
slope = trace[:slope]
intercept = trace[:intercept]
# Draw the line
plot([xmin, xmax], slope * [xmin, xmax] .+ intercept, color="black", alpha=0.5)
ax = gca()
ax[:set_xlim]((xmin, xmax))
ax[:set_ylim]((xmin, xmax))
end;
figure(figsize=(3,3))
render_trace(trace);
エクササイズ
同じアドレスを2回使うデータ生成関数を書けとのことなのでsloe * (x + x^2) + interceptの2次関数で試してみたいと思います。
@gen function line_model2(xs::Vector{Float64})
n = length(xs)
slope = @trace(normal(0, 1), :slope)
intercept = @trace(normal(0, 2), :intercept)
for (i, x) in enumerate(xs)
@trace(normal(slope * x + slope * x * x + intercept, 0.1), (:y, i))
end
return n
end
function render_trace2(trace; show_data=true)
xs = get_args(trace)[1]
xmin = minimum(xs)
xmax = maximum(xs)
if show_data
ys = [trace[(:y, i)] for i=1:length(xs)]
scatter(xs, ys, c="black")
end
slope = trace[:slope]
intercept = trace[:intercept]
xl = collect(range(xmin, stop=xmax, length=100))
yl = slope * xl .+ slope * map(x -> x^2, xl) .+ intercept
plot(xl, yl)
ax = gca()
ax[:set_xlim]((xmin, xmax))
ax[:set_ylim]((findmin(yl)[1], findmax(yl)[1]))
end
figure(figsize=(16,8))
traces = [Gen.simulate(line_model2, (xs,)) for _=1:12];
for (i, trace) in enumerate(traces)
subplot(3, 6, i)
render_trace2(trace)
end
2問目は周波数がgamma(5, 1)に従い、振幅はgamma(1, 1)の正弦波を書くものです。周波数、振幅はガンマ分布に従っておりGen.gmmma関数を使っています。ガンマ分布自体はトラフィックの待ち時間分布などに応用されているようでgamma(5, 1)は期間1の間に1度おきる事象が5回発生するまでにかかる時間の分布を表していて、実際にその時のガンマ分布を図示化すると以下のようになります。
import Pkg; Pkg.add("Distributions")
d=Gamma(5, 1)
xl = collect(range(0, stop=10, length=50))
plot(xl, map(x -> cdf(d, x), xl))
plot(xl, map(x -> pdf(d, x), xl))
グラフより4の辺りが確率のピークになっていることが確認できます。
あとは正弦波の位相のずれをphage、周波数をperiod、振幅をamplitudeとした時、座標xに対するy座標はamplitude * sin(2 * pi * x / period + phase)となるので、これまでと同様の手順で図示化できます。
今回は先にモデルを決めた状態でy座標を生成しましたが、次は実際に必要となる座標のデータからモデルを発見する手順をみていきたいと思います。
データからモデルを推論する
学習データとして以下を準備します。
xs = [-5., -4., -3., -.2, -1., 0., 1., 2., 3., 4., 5.]; ys = [6.75003, 6.1568, 4.26414, 1.84894, 3.09686, 1.94026, 1.36411, -0.83959, -0.976, -1.93363, -2.91303]; figure(figsize=(3,3)) scatter(xs, ys, color="black");
データからモデルを生成する場合は(Gen.importance_resampling)https://probcomp.github.io/Gen/dev/ref/inference/#Importance-Sampling-1を使います。この関数に対してデータ生成関数と入力、出力の教師データ、それから計算回数を渡します。以下の関数を準備しておけば
function do_inference(model, xs, ys, amount_of_computation)
# Create a choice map that maps model addresses (:y, i)
# to observed values ys[i]. We leave :slope and :intercept
# unconstrained, because we want them to be inferred.
observations = Gen.choicemap()
for (i, y) in enumerate(ys)
observations[(:y, i)] = y
end
# Call importance_resampling to obtain a likely trace consistent
# with our observations.
(trace, _) = Gen.importance_resampling(model, (xs,), observations, amount_of_computation);
return trace
end;
先に準備していた線形モデルのデータ生成関数と教師データ、計算回数を渡すことでモデルが生成できます。do_inference(line_model, xs, ys, 10)の最後の引数として渡している計算回数を変えることでモデルがよりフィットすることが確認できると思います。
trace = do_inference(line_model, xs, ys, 10) figure(figsize=(3,3)) render_trace(trace);
とりあえずデータ生成関数さえ準備できれば、あとは何度か計算してよかった結果を使えば良いのかと思いました。 線形モデルについても同様の手順で行うことができます。学習デーがが以下の場合
xs = [-5., -4., -3., -.2, -1., 0., 1., 2., 3., 4., 5.]; ys_sine = [2.89, 2.22, -0.612, -0.522, -2.65, -0.133, 2.70, 2.77, 0.425, -2.11, -2.76]; figure(figsize=(3, 3)); scatter(xs, ys_sine, color="black");
先に用意したsine_model関数を使って以下のように求められます。
trace = do_inference(sine_model, xs, ys_sine, 1000) figure(figsize=(3,3)) render_sine_trace(trace);
ここまで試してきた限りだと他の機械学習ライブラリが学習を繰り返してモデルの適合率をあげるアプローチをしているのに比べて、Genはパラメータの確率分布モデルを先に決めておいて何度か計算をしたうえで最適なもの選択するので過学習を気にする必要はないのかという気がしてその辺りが統計を知らない人向けでもあるのかと思いました。ただ先にデータ生成モデルを準備しておく必要があり、この辺りは確率的プログラミングを学んでおく必要がありそうで、その辺りがわからないとGenが初心者向けなのかわからない気がしました。
Pythonでの確率的プログラミングのライブラリではEdwardがあるらしいです。こちらは確率分布モデルに対してベイズ推論でアプローチするらしいです。 Juliaの方も確率分布モデルを推論する仕組みがありそうで、よく調べきれていませんが以下のtutorialではデータ生成関数の推論が行えるように見えました。
probcomp.github.io これまでの機械学習と比べて学習の対象が確率分布モデルになっているような気がします。
JuliaとかEdwardとか人の名前が多くて検索しずらくなるかも
