2018-03-22

【数理統計】ラオーブラックウェル化（Rao-Blackwellization）

統計

理解したいのはこっちでした。

本質的には同じ問題だけど。

ラオーブラックウェル化（Rao-Blackwellization）は、条件付き期待値を利用して期待値を計算する方法です。

サンプルが独立なら、ブラックウェルーラオの定理より直接期待値を計算するより精度が良くなる事が保証されます（不思議）。

名前の順番がひっくり返っているのは知りません（書籍でこう書いてあったからそのまま書いてます）

ブラックウェルーラオの定理の証明はこちら

hotoke-x.hatenablog.com

ラオーブラックウェル化（Rao-Blackwellization）

確率変数 $\boldsymbol{x}$ を二つのブロックに分割し、 $\boldsymbol{x} = (\boldsymbol{x_1}, \boldsymbol{x_2}) (x_i \in \chi_i)$ と表す。このとき、 $x_1$ の統計量の期待値

f:id:hotoke-X:20180322193210p:plain

を求めることを考えます。

$\pi (\boldsymbol{x_1}, \boldsymbol{x_2})$ からサンプリングされた $\left( \boldsymbol{x^{(1)}}, \ldots, \boldsymbol{x^{(t_0 + T)}} \right)$ を用いれば、上式の期待値は以下のように計算できます。

$\displaystyle \hat I = \frac{1}{T} \sum_{t=1}^{T} g(\boldsymbol{x_1^{(t_0 + T)}})$

ここで、 $\boldsymbol{x^{(t)}} = \left(\boldsymbol{x_1^{(t)}}, \boldsymbol{x_2^{(t)}} \right)$ であることに注意。

さらに $\boldsymbol{x_2}$ の条件付き分布で同時分布を分解すると

f:id:hotoke-X:20180322194353p:plain

となるので、この式の真ん中部分

$\displaystyle h(\boldsymbol{x_2}) = \boldsymbol{\mathrm{E}} \left[g(\boldsymbol{x_1}) | \boldsymbol{x_2} \right] = \int_{\chi_1} g(\boldsymbol{x_1})\pi(\boldsymbol{x_1} | \boldsymbol{x_2}) \mathrm{d}\boldsymbol{x_1}$

を解析的に求められれば、

$\displaystyle \hat I_{RB} = \frac{1}{T} \sum_{t=1}^{T} h(\boldsymbol{x_2}^{(t_0 + t)})$

によって計算できる。

なんだ面倒くさいと思うが、実はラオーブラックウェルの定理から

$\displaystyle \boldsymbol{\mathrm{V}} (\hat I) \geq \boldsymbol{\mathrm{V}} (\hat I_{RB})$

が成立し、より良い推定値になり得ることがわかる。

参考書籍

応用を目指す数理統計学（国友直人）
ベイズ計算統計学（古澄英男）

2018-03-22

【数理統計】ブラックウェルーラオ（Blackwell-Rao）の定理

統計

はい。

私がちゃんと理解したかった定理です。

復習したら見た瞬間理解したのですが、一応メモっておきます。

初学者向けの説明がなかなかネットに落ちていなかった（気がする）ので。

初めに行っておくと、この定理は

「十分統計量で条件づけた不偏推定量の分散は、他の不偏推定量の分散より大きくなることはない。」

って定理です（違ったらコメントください）。

十分統計量に依存した不偏推定量を使っとけばとりあえず無難ですよってことですな。

数学アレルギーでも、この程度の認識は持っておいた方が良いでしょう。

十分統計量についてはこちら hotoke-x.hatenablog.com

不偏推定量についてはこちら hotoke-x.hatenablog.com

証明

証明は結構簡単です。

今、未知母数の推定量は不変推定量しか考えないとします。このとき、標本 $\boldsymbol{X}$ について $\theta$ の不偏推定量 $\phi (\boldsymbol{X})$ 、十分統計量 $T(\boldsymbol{X})$ を考えます。この時十分統計量で条件付けした統計量を $\varphi (t) = \boldsymbol{\mathrm{E}} \left[\phi (\boldsymbol{X}) | T = t \right]$ とすると

$\boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \theta \right)^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right) + \left(\varphi (T) - \theta \right) \right]^2 \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right)^2 \right] + \boldsymbol{\mathrm{E}} \left[ \left(\varphi (T) - \theta \right)^2 \right] \\ \geq \boldsymbol{\mathrm{V}} \left[ \varphi (T) \right]$

十分統計量に依存しない不偏推定量を用いた場合、下から2段目の式の第一項の分、推定量が悪くなり得ることを示しています。

まとめるとブラックウェルーラオ（Blackwell-Rao）の定理は以下のようになります[1]。

標本 $\boldsymbol{X} = (X_1, \ldots, X_n)'$ の各要素が互いに独立同分布に従うとき（independent and identically distributed, i.i.d.と書くことも多い）、 $T(\boldsymbol{X})$ を十分統計量、統計量 $\phi (\boldsymbol{X})$ は母数 $\theta$ の不偏推定量とする。このとき

$\varphi (t) = \boldsymbol{\mathrm{E}} \left[ \phi (\boldsymbol{X}) | T=t \right]$ は不変推定量
$\boldsymbol{\mathrm{V}} \left[ \varphi (T) | \theta \right] \leq \boldsymbol{\mathrm{V}} \left[\phi (\boldsymbol{X}) \right]$ となる

参考書籍

応用を目指す数理統計学（国友直人）

2018-03-22

【数理統計】最尤推定量とバイアス

統計

尤度

尤度については

hotoke-x.hatenablog.com

を参照。

例によって離散か連続かはあまり気にせず書きます。

尤度関数とは、母集団分布からの独立標本が得られたとき、確率関数 $p(x|\theta)$ から同時確率関数を母数 $\theta$ の関数と見た

$$L_n (\theta | x_1, \ldots, x_n) = \prod_{i=1}^{n} p(x_i | \theta) $$

のことでした[1]。

この尤度関数を最大にするような $\theta$ の値を最尤推定値といい、その推定量 $\hat \theta$ を最尤推定量（maximum likelihood estimator, MLE）と呼びます。私の読むレベルの文献では $\hat \theta_{ML}$ と書かれていることが多いです。

最尤推定

最尤推定すれば未知母数の推定もうまくいきそうですが本当にそうでしょうか。百聞は一見に如かずということで確かめてみます。

ここで、正規分布の平均と分散を未知として、最尤推定する問題を考えます。

母集団分布を平均 $\mu$ 、分散 $\sigma^2$ の正規分布 $\mathcal{N}(\mu, \sigma^2)$ とする。標本を独立にn個サンプリングし、観測値 $x_1, \ldots, x_n$ が得られたとするとその尤度関数は

$\displaystyle L_n (\mu, \sigma | x_1, \ldots, x_n) = \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right)^n \exp \left\{ -\frac{1}{2\sigma^{2}} \sum_{i=1}^{n} (x_i - \mu)^2 \right\}$

で与えられ、その対数尤度関数は

$\displaystyle l_n (\mu, \sigma) = -\frac{n}{2} \log (2\pi) - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \sum_{i=1}^{n} (x_i - \mu)^2$

となります。

母数 $\mu$ で微分して、 $\mu$ について解けば母数 $\mu$ についての最尤推定量が得られます。

$\displaystyle \frac{\partial l_n}{\partial \mu} = \frac{1}{\sigma^2} \sum_{i=1}^{n} (x_i - \mu) = 0$

$\displaystyle \hat \mu_{ML} = \frac{1}{n} \sum_{i=1}^{n} x_i = \bar x$

この推定量を使って、分散も推定ができます。そしてこの推定量を使って他の推定量を得るという作業に落とし穴があるので注意（後述）。

対数尤度に $\mu$ の最尤推定量 $\hat \mu_{ML}$ を代入して母数 $\sigma$ で微分すると

$\displaystyle \frac{\partial l_n}{\partial \sigma} = - \frac{n}{2\sigma^2} + \frac{\sum_{i=4}^{n} (x_i - \hat \mu_{ML} )^2}{2\sigma^4} =0$

が得られ、 $\sigma$ について解くと

$\displaystyle \hat \sigma_{ML}^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar x)^2$

が得られます。

以上の手順で、正規分布の母数（平均と分散）の最尤推定量

$\displaystyle \hat \mu_{ML} = \frac{1}{n} \sum_{i=1}^{n} x_i = \bar x$

$\displaystyle \hat \sigma_{ML}^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar x)^2$

が得られました。

バイアス

先ほどの最尤推定量を得る過程で推定量を使って他の推定量を得るなんて少し怪しいことをしたわけですが、怪しいという直感は当たっていて、分散は不偏推定量になっていません。不偏というのは推定量が真の母数の周りに対称に分布することを言います。つまり、不偏でない推定量は母数を過小（あるいは過大）評価することを意味します。

推定量と母数の二乗誤差の期待値（mean square error, MSE）を計算して見ると

$\boldsymbol{\mathrm{E}} \left[ |\hat \theta - \theta|^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( (\hat \theta - \boldsymbol{\mathrm{E}} (\hat \theta)) + (\boldsymbol{\mathrm{E}} (\hat \theta) - \theta) \right)^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( \hat \theta - \boldsymbol{\mathrm{E}} (\hat \theta) \right)^2 \right] + \left[ \left( \boldsymbol{\mathrm{E}} (\hat \theta) - \theta \right)^2 \right] \\= \boldsymbol{\mathrm{V}} \left( \hat \theta \right) + \left( \boldsymbol{\mathrm{E}} (\hat \theta) - \theta \right)^2$

となり、右辺第二項の分推定量が偏ることがわかります。不偏推定量であればこのバイアス項は0になります。数学アレルギー勢にために一応説明しておくと、 $\boldsymbol{\mathrm{V}}$ は分散（Variance）という意味です。

先ほどの正規分布の例で、分散の最尤推定量が不偏推定量になっていないというのはこのバイアス項が0にならないことを意味しています。

それは困ったと思うかもしれませんが、どの程度偏っているかわかっているなら補正すればよいだけです。

分散の最尤推定量の期待値を計算してみると（計算は気が向くか、コメントいただければまた書きます（mathjax疲れた））

$\displaystyle \boldsymbol{\mathrm{E}} (\hat \sigma_{ML}^2) = \frac{(n-1)}{n} \sigma^2$

となり、 $(n-1)/n$ 倍だけ過小評価されている事がわかります。

じゃあ最尤推定値をこの分だけ補正して母数を推定すれば良いのです。

統計で習う分散の式が

$\displaystyle \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \mu)^2$

となっているのはこの為です。

式だけでいきなり「ほぉそうか」と納得できる人はこれで良いのですが、数学アレルギー勢は釈然としないと思います。

絵で描くとわかりやすいので、気が向いたら絵でも描こうと思います。

今日は疲れたのでこの辺で。

参考書籍

応用を目指す数理統計学（国友直人）

2018-03-22

【数理統計】尤度関数と十分統計量

統計

勉強したことを自分なりにまとめたものです。
間違い等があれば指摘していただけると嬉しいです。

私は「せいきぶんぷ～？」、「のんぱらめとりっく～？」という状態から

応用を目指す数理統計学（国友直人）

で勉強しました。初学者が手を出すには敷居が高いですが、（本気で）頑張ればなんとかなるように書かれています。参考文献も充実しているのでとりあえず持っておくのはおすすめです。安いし。（アフィリエイト等はないので安心して買ってください。）

以下、応用を目指す数理統計学（国友直人）に基づいています。

自分なりの解釈も書いたりしてます。

Chapter 7 統計的推論より

統計量・推定量・推定値

統計量とは「標本の関数」のこと
推定量とは「統計量を使って母集団を表現している未知母数を推測するときの統計量のこと」
推定値とは「実際にデータを代入した場合の推定量のこと」

標本の関数でさえあれば統計量と呼ぶようです。その中でも、母数を推測する統計量を推定量と区別することに注意が必要といったところでしょうか。

「応用を目指す数理統計学」では正規分布に従うことの妥当性を標本積率を使って調べる方法を例に挙げてくれています（p108, 109）。（書きすぎるとただの本のコピペになるので、気になる方は購入してください）。母集団と標本の積率を対応づけるので積率法と呼ばれるものの一例のようです。最近では積率法を一般化した一般化積率法がよく使われているとのこと。10章で解説してくれているそうなのでまあたとで勉強しよう。

尤度関数と十分統計量

尤度関数とは「母数の関数」のこと
十分統計量とは「母数に関する情報損失が無い統計量」のこと

尤度関数

簡単のため、離散確率分布の場合を考えます（連続でも同様）。標本は大文字、確率変数は小文字で書くのがスタンダードっぽい。

$\theta$ を母数、 $x$ を標本の確率変数、 $p(x|\theta)$ を母集団が従う確率変数とする。

この時、 $n$ 個の標本 $X_1, X_2, \ldots , X_n$ の同時確率関数は

$$ p(x_1, x_2, \ldots , x_n | \theta) = \prod_{i=1}^{n} p(x_i|\theta) $$

で与えられる。これを、同時確率分布ではなく、母数 $\theta$ の関数として考えると

$$L_n (\theta | x_1, \ldots, x_n) = \prod_{i=1}^{n} p(x_i | \theta) $$

と書ける。これを尤度関数と呼ぶ。 $L$ はLikelihoodの頭文字。

得られた標本から計算された尤度関数の意味としては「起きた事象はどの程度起こりやすい事象だったのか」を定量化したことになりそうです。俗にいう「尤もらしさ」とはこのことらしい。ただし、既知にしろ未知にしろ、「ある分布に関する母数」を設定してあるということはパラメトリックな分布を想定していることになります。真の分布ではなく、統計モデルに落とし込んでいるということは意識として持っておく必要がありますね。

実際には、（生物系では）次々にデータ（標本）が増えるシチュエーションはなかなかないので、標本は固定で $\theta$ を色々動かして母数を推定することになります（例えば最尤推定）。

十分統計量

データから母数 $\theta$ を推定する事をパラメトリック推測と呼ぶ。この時使う統計量が、母数に関して情報損失が無い場合、その統計量を十分統計量という。すなわち、十分統計量は母数 $\theta$ とは無関係の関数となる。

Fisherの因子分解定理が成立すれば十分統計量ってことらしい。

Fisherの因子分解定理
標本 $X_1, \ldots, X_n$ に対して同時密度関数 $f(x_1, \ldots, x_n)$ が存在し、 $\boldsymbol{x} = (x_1, \ldots, x_n)$ とする。この時

$$f (\boldsymbol{x} | \theta) = g(T(\boldsymbol{x}) | \theta) h(\boldsymbol{x})$$

と分解できれば $T$ は十分統計量である。ということらしい。

証明は非常に簡単で、 $f(\boldsymbol{x} | \theta)$ にこの分解定理を適用して確率密度を計算しなおしてみればよい。

$$f (\boldsymbol{x | \theta}) = f (\boldsymbol{x} | T=t)g(t|\theta) = h(\boldsymbol{x}; t)g(t|\theta)$$

ここでセミコロンは、変数かgivenな値か区別する記号（だと思われる）。これを全事象について積分して割ってやれば確率密度を計算できるので、全事象を $A_t = \{\boldsymbol{x} | T(\boldsymbol{x}) = t \}$ として

$$\int_{A_t} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} = g(T(\boldsymbol{x}) | \theta) \int_{A_t} h(\boldsymbol{x}; t) \mathrm{d}\boldsymbol{x}$$

とすると、

$$f(\boldsymbol{x} | T=t(\boldsymbol{x})) = \frac{f (\boldsymbol{x | \theta}) }{\int_{A_t} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} } = \frac{h(\boldsymbol{x}; t)g(t|\theta)}{g(t | \theta) \int_{A_t} h(\boldsymbol{x} ; t) \mathrm{d}\boldsymbol{x}} = \frac{h(\boldsymbol{x}; t)}{\int_{A_t} h(\boldsymbol{x} ; t) \mathrm{d}\boldsymbol{x}}$$

となって母数 $\theta$ と無関係になる。

言われてみれば当然なんだけど、十分統計量だけ調べれば母数の推測は不要という点で強力ですね。

2018-03-08

Julia tips #9: CVODE (Sundials.jl) の数値積分が終わらない時の回避方法

Julia

追記
どうやら、せっかく設定したhminは内部でフラグチェックが行われないようです。
solve関数のmaxitersオプションで対応するしかないかもしれません。また後程検証してみます。

以下は検証が十分ではありません。
ソルバ自体の動作に支障はありませんが、期待通りに動いているかは検証する必要があります。
そのことをご理解いただいたうえでご覧ください。

2018/03/08当時の環境
- Julia 0.6.2 - Sundials 1.2.0

事の発端

ODEモデルのパラメータ推定をしていたら、非線形なモデルだとどうしても途中でソルバがフリーズしてしまう...。

原因はソルバが無限に数値積分のステップを続けてしまうことでした。

数値計算に疎いのでどうしてそんなことになってしまうのかはちゃんとわかっていないのですが、stiffなODEを解こうとすると、数値積分時に細かいステップを刻まないといけなくなるのが原因です。

「DifferentialEquationsから簡単にSundialsが使える！」と喜んだのもつかの間。

DifferentialEquations.jlでもこの件について触れられていました。

Frequently Asked Questions · DifferentialEquations.jl

公式によると、現状アルゴリズムを変えたり、ソルバを変えたりするしかないようです。

しかし、stiffなODEが解けないような状況においてCVODE以外ではお話にならないのです（経験上）。

そもそもODEが解けないようなパラメータは捨ててしまえばよいので、途中で計算を止めてなんとか捨てたいわけです。

Sundials? CVODE?

Sundialsは常微分方程式のソルバ、CVODEを作っているところです（他に幾つも作ってる）。多用しないのでわかりませんが、良く名前を聞くのでシェアは大きいと思います、MATLABにもあります。

stiffな問題をCVODEに解かせたら他のソルバより圧倒的に高速で数値積分してくれます。
stiffでなくてもめっちゃ速いです。過去記事中のリンク先でJuliaのソルバで速度比較をした方がいます。

hotoke-x.hatenablog.com

ちなみにCVODEとよく比較されるソルバにLSODAっていうのがあります。自分の場合はCVODEの方が5倍くらい速かったです。問題にもよると思うので、CVODEがダメだったら試す価値はあるかもしれません。JuliaならPkg.add("LSODA")ってやればDiffEq.jlのインターフェースから使えるようになります。以下のブログ記事ではLSODAも良いぞ！って言ってますね。

Solving Stiff ODEs

CVODEが途中で止まる

最強のCVODEといえど限界はあります。積分のステップが小さすぎると計算が終わらなくなります。途中でエラー吐けよって思うんですが、デフォルト設定だと無限にステップを刻めることになっています。こいつをちょちょいといじれば良いのです。

以下のhminというオプションをいじる（CVODE公式ドキュメント, p37より）↓ f:id:hotoke-X:20180308031520p:plain

The default value is 0.0となっていますね。

これがいかんわけです。

Sundials.jlを書き換える（自己責任で）

DiffEq.jlのインターフェースには残念なことに、hminを設定するオプションがありません。しかし、Sundials.jl自体ではちゃんとラップしてくれているのです。じゃあ直接Sundials.jlを使えばいいんですが、exampleを見てるとわかるんですがjuliaの可読性が損なわれます。何より面倒です。

なので、DiffEq.jlのインターフェースから使えるようにSundials.jlを少しだけいじりました。

自分でやる場合は自己責任でお願いします。 ...といっても、ミスっても書き換えたところを戻せば元通りです（2敗）

DiffEq.jlのインターフェースとなるsolve.jlをSundialsのパッケージディレクトリから探します。 f:id:hotoke-X:20180308033027p:plain

hminを引数としてsolveに追加（後ろにkwargs...があるので多分これやらなくても動くけど）。

f:id:hotoke-X:20180308033609p:plain

以下のようにhminをCVODEの設定として反映させる（ハイライト部分）。

f:id:hotoke-X:20180308033746p:plain

これで完了。

SundialsはCのライブラリなんですが、Sundilas.jlの作者さんがきれいにラップしてくれているおかげでインターフェース部分を少し書き換えるだけで良いです。作者さんありがとうございます。

後は動作確認したらよいです。

パッケージを書き換えたので、初めにプリコンパイルが走ります。

（以下ではCVODEが固まるほどのstiffなODEがどんなものかわからないので、ソルバが動くかどうかだけ確かめています。）

using DifferentialEquations
using Sundials
using Plots

function parameterized_lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
alg = CVODE_BDF()
prob = ODEProblem(parameterized_lorenz, u0, tspan, p)
sol = solve(prob, hmin=1.0e-8)

plot(sol,vars=(1,2,3))

以下のようにグラフがプロットされればとりあえず正常に動作しています。

f:id:hotoke-X:20180308043930p:plain

しかし、本当に設定が反映されたのだろうか...

（どなたか、CVODEをフリーズさせる程のめちゃめちゃ硬いODEを教えてください）

2018-02-17

Julia tips #8: パッケージの追加に失敗したとき

Julia

自宅PCにDifferentialEquations.jlを追加しようとしたら以下のようなエラーが...

julia> Pkg.add("DifferentialEquations")
ERROR: resolve is unable to satisfy package requirements.
  The problem was detected when trying to find a feasible version
  for package DiffEqPhysics.
  However, this only means that package DiffEqPhysics is involved in an
  unsatisfiable or difficult dependency relation, and the rootof
  the problem may be elsewhere.
  (you may try increasing the value of the JULIA_PKGRESOLVE_ACCURACY
environment variable)

パッケージのアップデートしたら解決するだろうと

Pkg.update()

もしてみたが解説せず。

同じことで困っている人を探したら以下のような投稿を発見。

discourse.julialang.org

いっぱいやり取りしてくれてますが、要は

julia>ENV["JULIA_PKGRESOLVE_ACCURACY"] = 10
julia>Pkg.resolve()

とするだけ。

何が起きてるのかはよくわかりませんが、バージョンの依存関係をいい感じに直してくれてるみたい？

後はPkg.addで好きなパッケージを追加していきましょう。

2018-02-12

Julia tips #7: Juliaでプロットのlegendを手動で変更したい

Julia

Juliaのプロットにて、legendの変更で躓いたので色んなバックエンドで簡単に試してみた。

Plotlyでのプロットが消えてしまいましたが、まぁいいでしょう。

~~legendをプロットの外へ移動したい場合はどうすればいいんだろうか...~~

解決しました。

これではだめですか？ https://t.co/KTMiog5fvY
— ktrst (@ktrst) 2018年2月12日

@ktrstさん、ありがとうございます。

2018-02-08

Jupyter notebookのフォントを変える（Windows 10）

Juliaで下付き文字を使うときにあまりにも見にくかったので、フォントを変えて多少なり誤魔化したいというのが事の発端。

幾らでも情報は転がっていますが、ぱっと調べた感じだとWindowsでいじっているのを見かけなかったので、念のためのメモ書きです。

以下のように記述したファイルをcustom.cssという名前で.jupyter/custom以下に保存します。

.jupyterファイルは大抵の人はユーザフォルダ直下にあると思います。

.CodeMirror pre {
    font-family: Consolas;
    font-size: 12pt;
}

WindowsなのでConsolasフォントを指定しています。

好きなフォントにしたらいいと思います。

また、フォントサイズを少し大きくして多少誤魔化しています。

2018-02-01

Windows 10でJulia専用のPython仮想環境を作る

Julia

Juliaでは、Pythonを呼び出して使うことができます。しかし、何も考えずに使おうとするとMinicondaのPython 2.7がインストールされてしまいます（2018/02/01現在）。

そこで、Anacondaで作ったJulia専用の仮想環境を作って、その中でJulia用のPythonパッケージ管理をしようという話です。

以下のような人を対象としています。

既にPythonをインストールしてあるのでそれを使いたい
既にあるPython環境を汚したくないのでJulia用の仮想環境を作って使いたい

Anacondaをダウンロードしてインストール

既にAnacondaをインストールしている人は読み飛ばしてください。

以下からWindows用のAnacondaインストーラをダウンロードしてインストールします。 https://www.anaconda.com/download/

インストールの仕方は公式のこちらで確認できます。 https://docs.anaconda.com/anaconda/install/windows

インストールの際、AnacondaのPythonをPathへ追加するかどうかのチェック画面が現れます。

WindowsにはデフォルトでPythonが入っていないので追加しちゃってよいでしょう。

f:id:hotoke-X:20180201203308p:plain 両方にチェックを入れときましょう。
（嫌な人はAnacondaを起動するときは、スタートメニューからAnaconda プロンプトを起動してください。）

Julia専用のAnaconda仮想環境の構築

コマンドプロンプト（またはAnacondaプロンプト）にて、以下のように入力して仮想環境を作ります。
ついでにjupyterも一緒にインストールしておきましょう。

C:\Users\username>conda create -n env_name python=3 jupyter

env_nameに自分の好きな仮想環境名を入れてください。

次に、activateして仮想環境の中に入ります。

また、Juliaで使いそうなパッケージをインストールしておきます。

C:\Users\username>activate env_name
C:\Users\username>conda install numpy pandas scipy matplotlib

念のため、Jupyter notebookのカーネルに作った仮想環境を追加しておく。

（多分やらなくても良い）

※必ず仮想環境に入った状態で行ってください。

C:\Users\username>ipython kernel install --user --name=env_name --display-name=env_name

仮想環境のPythonのパスを確認

C:\Users\username>where python
C:\Users\username\Anaconda3\envs\env_name\python.exe
C:\Users\username\Anaconda3\python.exe

仮想環境とAnacondaインストール時に入るPythonのパスが表示されます。

envsの下にあるpythonが仮想環境のpythonなのでパスをコピーしておきます。　　

仮想環境から抜けて、jupyterのパスを取得します。

C:\Users\username>deactivate
C:\Users\username>where jupyter
C:\Users\username\Anaconda3\Scripts\jupyter.exe

表示されたパスを環境変数JUPYTERへ追加します。

コントロールパネル＞システムとセキュリティ＞システム＞システムの詳細設定＞環境変数(N)

で環境変数の設定画面に行けます。

上段のユーザー環境変数に新規でJUPYTERという名前で環境変数を作成し、jupyterのパスを入力します。

Juliaが勝手に環境変数を読んでjupyterを見つけてくれます。

jupyter.exeまで入力しないと無視されるので注意

以上でPython側の準備は終わりです。

Juliaをインストール

公式からWindows用のJuliaをダウンロードしてインストール。

インストール場所はどこでも良いです。

https://julialang.org/downloads/

このままだとPCがjuliaを見つけられず

C:\Users\username>julia

と入力しても何も起きません。パスを教えてやる必要があります。

インストールしたフォルダ/binをユーザー環境変数Pathへ追加。

（julia.exeがあるのでわかりやすいと思います）

juliaのパッケージをインストールするフォルダの場所をデフォルトから変えたい場合は、環境変数JULIA_PKGDIRを作って好きなパスを入力して下さい。

追加の仕方は環境変数JUPYTERの設定の時と同じです。

次に、コマンドプロンプトを立ち上げなおしてjuliaを起動します。

まず初めにJuliaの環境変数を書き換えて、先ほど作った仮想環境のPythonの場所を教えておきます。

julia>ENV["PYTHON"] = "C:/Users/userneme/Anaconda3/envs/env_name/python.exe"

juliaの設定を初期化します。

PyPlotを入れるとPyCallが一緒についてきて、その際に仮想環境のpythonを見つけて使ってくれるようになります。

julia>Pkg.init()
julia>Pkg.add("IJulia")
julia>Pkg.add("PyPlot")
julia>using PyCall

以上で仮想環境のpythonを使うことができました。

まとめ

環境が汚れなくてハッピー！
Julia自体でも仮想環境を作れたら楽。

2018-01-31

Julia tips #6: Multiple-Try Metropolisを実装してみた

Julia 統計

Multiple-Try Metropolis (MTM)って？

メトロポリス・ヘイスティングスの多数決バージョンです。
今いる場所から次の候補を複数生成し、その値周辺が今いる場所より尤もらしければ候補の中の1つへ移動します。

細かい話を書く元気がないので、理論的な背景は気が向いたらまたまとめます。
というわけで実装して行きましょう。

MTMアルゴリズム

任意の提案関数 $T(x,y)$ を定義し $T(x, y) > 0 \Longleftrightarrow T(y, x) > 0$ を満たすとする。
また、 $\lambda (x,y)$ を $\lambda (x,y) > 0$ で対称な任意の関数であるとする。
マルコフ連鎖の現在の状態を $\boldsymbol{x_t}$ とおくと、MTMアルゴリズムは次のようになる。

$T(\boldsymbol{x_t}, \cdot)$ から $\boldsymbol{y_1}, \dots, \boldsymbol{y_k}$ を独立に生成する
重み $w(\boldsymbol{y_j}, \boldsymbol{x_t}) = \pi(\boldsymbol{y_j})T(\boldsymbol{y_j}, \boldsymbol{x_t})\lambda(\boldsymbol{y_j}, \boldsymbol{x_t})$ に従って ${\boldsymbol{y_1}, \dots, \boldsymbol{y_k} }$ から $\boldsymbol{y}$ をサンプリングして計算する。
$T (\boldsymbol{y}, \cdot)$ から $\boldsymbol{x_1^*}, \dots, \boldsymbol{x_{k-1}^*}$ をサンプリングし、 $\boldsymbol{x_k^*}=\boldsymbol{x_t}$ とする。
$y$ を以下の確率で受容する。 $\alpha = \mathrm{min} \left\{1, \frac{\Sigma_{k} w(\boldsymbol{y_k}, \boldsymbol{x_t})}{\Sigma_{k} w(\boldsymbol{x_k}, \boldsymbol{y})} \right\}$

今回はベイズ計算統計学（統計解析スタンダード）p93 例3.7を実装してみた。提案関数と対称関数は以下のようになっている。

$$T(x, y) = x_t + \mathcal{N}\left(0, 5\sqrt{2}\right)$$ $$\lambda(\boldsymbol{y}, \boldsymbol{x})=\left(\frac{T(x, y)+T(y, x)}{2}\right)^{-1}$$

目的関数は、

$$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$$

である。

使うパッケージのインポート

using Distributions
using StatsBase
import Plots; plt = Plots; plt.pyplot()
using LaTeXStrings

# グラフ周りの設定
fntsm = plt.font("serif")
fntlg = plt.font("serif", 20)
plt.default(titlefont=fntlg, guidefont=fntsm, tickfont=fntsm, legendfont=fntsm)
plt.default(size=(800, 600))

目標分布の定義

$$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$$

const μ = [-5.0, 5.0]
const σ = [1.0, 1.0] 

function target(x::Array{Float64, 1})
    arr = zeros(length(x))
    Normal_val_1(x) = 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    for i in 1:length(x)
        arr[i] = 1/3 * Normal_val_1(x[i]) + 2/3 * Normal_val_2(x[i])
    end
    
    return arr
end

function target(x::Float64)
    Normal_val_1(x) = 0.45 * 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 0.45 * 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    val = 1/3 * Normal_val_1(x) + 2/3 * Normal_val_2(x)
    return val
end

x = Array{Float64}(linspace(-10, 10, 1000))
y = target(x)
l = plt.@layout [a b; c{0.2h}]
eq = L"$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$"
p = plt.plot(x, y, title="Target distribution", lab=eq)

f:id:hotoke-X:20180131003314p:plain

実装

まず、MTMに必要な関数やら値を保持する型を作っとく。

mutable struct Sampler
    x_present::Float64   # 現在の値
    target::Function   # 目標分布関数
    propose_next::Function   # 提案関数
    q::Function   # (x_t+1, x_t)の順で渡す   # 推移核（上で書いたT(x, y)）
    λ::Function   # 対称関数
    burn_in::Int   # 焼きなましの回数
    count::Int   # 何回目の更新か保持
    random_state::MersenneTwister   # 乱数ジェネレータ
end

MTM本体

module MTM
using StatsBase   
    function present(sampler)   # 現在の値を返す関数
        return sampler.x_present
    end

    function next(sampler, n_candidate)   # MTM
        w_x2y = zeros(n_candidate)
        w_y2x = zeros(n_candidate)
        proposed = zeros(n_candidate)
        proposed_r = zeros(n_candidate)
        accept = 0
    
        # Σⱼ w(yⱼ,xₜ)
        for i in 1:n_candidate
            # q(yⱼ|xₜ)でyⱼをサンプリング
            proposed[i] = sampler.propose_next(sampler.x_present)
        
            # q(yⱼ|xₜ) と q(xₜ|yⱼ) を計算
            q_x2y = sampler.q(proposed[i], sampler.x_present)
            q_y2x = sampler.q(sampler.x_present, proposed[i])
        
            # w(yⱼ, x) = π(yⱼ)q(yⱼ|xₜ)λ(yⱼ, x)
            w_x2y[i] = sampler.target(proposed[i]) * q_x2y * sampler.λ(q_x2y, q_y2x)
        end
        
    
    
        # w(y,xₜ)したがってyをサンプリング
        y = sample(proposed, Weights(w_x2y))
    
        # Σⱼ w(xₜ, yⱼ)
        for i in 1:(n_candidate-1)
            # q(xⱼ|y)でxⱼをサンプリング
            proposed_r[i] = sampler.propose_next(y)
        
            # q(xⱼ|y) と q(y|xⱼ) を計算
            q_y2x = sampler.q(proposed_r[i], y)
            q_x2y = sampler.q(y, proposed_r[i])
        
            # w(xⱼ, y) = π(xⱼ)q(xⱼ|y)λ(xⱼ, y)
            w_y2x[i] = sampler.target(proposed_r[i]) * q_y2x * sampler.λ(q_y2x, q_x2y)
        end
    
        # set xₘ=x
        proposed_r[end] = sampler.x_present
    
        # q(xₘ|y) と q(y|xₘ) を計算
        q_y2x = sampler.q(proposed_r[end], y)
        q_x2y = sampler.q(y, proposed_r[end])
        w_y2x[end] = sampler.target(proposed_r[end]) * q_y2x * sampler.λ(q_y2x, q_x2y) 
        
    
        # αの計算
        odds = sum(w_x2y)/sum(w_y2x)
        
        # 受容するか棄却するか決定
        if rand(sampler.random_state) < odds
            sampler.x_present = y
            accept = 1
            return y, accept
        else
            return sampler.x_present, accept
        end
    end

    function burn_in(sampler, n_candidate)   # 焼きなまし期間
        arr = zeros(sampler.burn_in)
        accept = zeros(sampler.burn_in)
        for i in 1:sampler.burn_in
            arr[i], accept[i] = next(sampler, n_candidate)
            sampler.count += 1
        end
        return arr
    end

    function run(sampler, n_steps, n_candidate)   #好きなだけ回す
        arr = zeros(n_steps)
        accept = zeros(n_steps)
        for i in 1:n_steps
            arr[i], accept[i] = next(sampler, n_candidate)
            sampler.count += 1
        end
        return arr, accept
    end

end

上で作った MTMを動かすところ。

const σ_proposal = 5*sqrt(2)  


const seed = 0
rng = srand(seed)   # 乱数シードの固定とジェネレータの取得
gaussian(xₜ) = Normal(xₜ, σ_proposal^2)   # 推移核の分布関数
generate(xₜ) = rand(gaussian(xₜ))   #  推移核の分布関数からサンプリングする関数
uniform = Uniform(-10, 10)   # 一様分布関数
inits = rand(uniform)   # 初期値の生成

target = target   # 意味ないけど、見た目のため
q(y, xₜ) =  pdf(gaussian(xₜ), y)
λ(x, y) = ((x+y)/2)^(-1) # MTM-inv scheme（下の参考文献参照）

const burn_in = 1000000
const num_sample = 100000
cnt = 0

sampler_MTM = Sampler(inits, target, generate, q, λ, burn_in, cnt, rng)

MTMの実行

ベイズ計算統計学（統計解析スタンダード）p93 例3.7の例と同じように、候補数を1, 5, 10と変えてみた。
候補数1は単純なメトロポリス・ヘイスティングス法になります。

焼きなまし

const n_candidate=10   # MTMでつかう候補の数（ここでは10）
samples_burn_in = MTM.burn_in(sampler_MTM, n_candidate)

本番。ステップごとのサンプルと受容フラグが返ってくる。

samples, accept = MTM.run(sampler_MTM, num_sample, n_candidate)

結果の確認

n_slice = 3000
lags =  collect(1:200)
auto_correlation = autocor(samples[end-n_slice:end], lags)


layout = plt.@layout [a b; c]

p1 = plt.bar(auto_correlation, title="Auto correlation",
                lab="", xlab="Lag", ylab="Auto\ncorrelation (A.U.)", 
                xlim=(-10,200), ylim=(-0.1,1))
p2 = plt.plot(samples[end-n_slice:end], title="Trajectory", lab="")
p3 = plt.histogram(samples , title="Target distribution", lab="", norm=true, bins=128)
p3 = plt.plot!(x,y, color="magenta", lab="")

p=plt.plot(p1, p2, p3, layout=layout)
print(sum(accept[end-n_slice:end]) / length(accept[end-n_slice:end]))   # 採択率を算出

f:id:hotoke-X:20180131010445p:plain — 候補数1（採択率：0.04798400533155615）

f:id:hotoke-X:20180131010726p:plain — 候補数5（採択率：0.20626457847384205）

f:id:hotoke-X:20180131010851p:plain — 候補数10（採択率：0.31756081306231254）

できました。
候補を増やすと、目標分布の2つの山を頻繁に行き来しているのがわかります。
採択確率も候補を増やしたほうが良くなっています。
ただ、候補の数分計算量が増えるので注意ですね。

まとめ

実装は簡単だったが、はてなブログに書くのが非常に疲れt
間違いがあったらご指摘ください。
ちなみに以前書いたメトロポリス・ヘイスティングスは少し間違ってます（笑）

参考文献

Pandolfi, S., Bartolucci, F. & Friel, N. A generalized multiple-try version of the Reversible Jump algorithm. Comput. Stat. Data Anal. 72, 298–314 (2014).