2020-12-23

Multioutput Gaussian Processで遊ぶ（Fun Advent Calendar 2020 #23）

機械学習

0．はじめに

Fun Advent Calender 2020の23日目です*1。僕はこのAdvent Calendarが作成された大学に所属する大学院生で、現在は博士後期課程3年（D3）に在籍しています。市内にある理学療法士・作業療法士養成校の教員（常勤）も務めており、研究と講義に追われる日々を送っています。研究では片麻痺の患者さんを対象として様々なデータをごにょごにょしています。専門分野はありません。
（※だいぶ前にTwitterの名前を2歳に変えてしまいましたがseibisiです。ちなみに2歳ではなく27歳です）

今回はMulti-output Gaussian Process（MOGP）というもので遊んでみました。MOGPはガウス過程の多出力版であり、正直あまり上手い日本語訳が分からない*2ので、英語のままにしておきます。Gaussian Processはそのままガウス過程であり、機械学習に興味がある人であれば1度は聞いたことのあるものだと思います。

この記事は、ガウス過程……聞いたことはあるけどなんだかよく分からない……難しい……という学部生を対象として、ガウス過程にHello Worldしてみよう！という目的で書いてみました。機械学習つよつよの方にとってはつまらない内容かもしれませんが、この記事を書くに当たって参考にした文献も載せておくので、理論に興味がある方はそちらを参照してください。

1．前準備

ガウス過程の話をする前に多変量正規分布について述べたいと思います。多変量正規分布はガウス過程における大きな立役者*3となります。この記事では終始一貫して多変量正規分布の話しか出てきません。

まずは定義から述べていきます。確率変数ベクトル $\displaystyle {\bf{x}}$ が平均ベクトル $\displaystyle {\boldsymbol{\mu}}$ 、分散共分散行列 $\displaystyle {\boldsymbol{\Sigma}}$ の多変量正規分布に従う時、多変量正規分布は次の式で定義されます。

$\displaystyle N( {\bf{x}} | {\boldsymbol{\mu}}, {\boldsymbol{\Sigma}} ) = \frac{1}{ (2\pi)^D \sqrt{ | {\boldsymbol{\Sigma}} | } } \exp ( -\frac{1}{2} ( {\bf{x}} - {\boldsymbol{\mu}} )^T {\boldsymbol{\Sigma}}^{-1} ( {\bf{x}} - {\boldsymbol{\mu}} ) )$

確率変数が複数存在する場合は確率変数間の対応関係を考える必要があり、この情報は分散共分散行列に存在しています。対応関係については、本学でも統計もしくは統計に関わる講義で共分散という言葉は1度は聞いたことがあるのではないでしょうか。聞いたことが無い人は、共分散と同じ意味を持つ言葉で相関という言葉が存在しているので、イメージとしてはそれでも大丈夫です。分散共分散行列は非常に重要なので、この記事を読む間は頭の片隅に置いておくと良いかもしれません。

2．Single-output Gaussian Process (SOGP)

まずは出力が1つのガウス過程（SOGP）から始めていく事にします。

2-1　SOGPの雰囲気

まずは下の図を見てみます。

図は何やら色とりどりのぐにゃぐにゃがいっぱい並んでいる様子です。横軸は $\displaystyle x$ 、縦軸には $\displaystyle f(x)$ と書いてあります。このぐにゃぐにゃたちの正体は、ガウス過程から生成されたサンプルです。え？グラフを生成？と混乱する人もいるかもしれませんが、次の図を見てみましょう。

図は、5、10、50、100変量の正規分布から乱数を発生させて結果を可視化したものです。横軸は確率変数ベクトル内の要素の番号（確率変数の番号）、縦軸はその実現値となるように可視化しました。なんとなくグラフっぽく見えるような気がします。確かに生成されているような…？

上の図は変量数を5、10、50…と進めていきましたが、これをさらに100、500、1000……さらには変量数が無限の正規分布からデータをサンプルして、同じようにグラフを描いたらどうなるでしょうか？これが1番最初の図になります。1つ1つのグラフは多変量正規分布からの「サンプル」と考える事もできると思います。

2-2　SOGPの定義

なんとなく雰囲気を掴めたところで、SOGPの定義を確認しておくことにします。出力をひとまとめにしたベクトルを $\displaystyle {\bf{f}} = ( f(x_1), f(x_2), ... ,f(x_N) )^T$ とした時、ガウス過程は次の式で定義されます。

$\displaystyle {\bf{f}}～GP( {\boldsymbol{\mu}}, {\bf{K}} )$

ここで、 $\displaystyle {\boldsymbol{\mu}}$ は平均ベクトル、 $\displaystyle {\bf{K}}$ は分散共分散行列であり、出力のベクトルは多変量正規分布に従います。話を簡単にするために、平均ベクトルを $\displaystyle {\bf{0}}$ としておきましょう。

$\displaystyle {\bf{f}}～GP( {\bf{0}}, {\bf{K}} )$

2-3　SOGPの分散共分散行列

さて、ガウス過程では分散共分散行列が大事（主観）という話をしたわけですが、SOGPの分散共分散行列についても見ていきましょう。SOGPでは分散共分散行列の $\displaystyle (i, j)$ 成分は $\displaystyle k(x_i, x_j)$ のように表現され、 $\displaystyle k(x_i, x_j)$ はカーネルと呼ばれています*4。カーネルには以下で述べるもの以外にも色々な種類があります。

ガウシアンカーネル

$\displaystyle k(x_i, x_j) = \theta_1 \exp ( -\frac{ |x_i - x_j| ^2 }{ \theta_2 } )$

指数カーネル

$\displaystyle k(x_i, x_j) = \exp ( -\frac{ |x_i - x_j| }{ \theta } )$

カーネルが何であるかの説明は詳しくはしませんが、SOGPでは分散や共分散を計算するための謎の式くらいでひとまず大丈夫です。試しにガウシアンカーネルと指数カーネルから生成したぐにゃぐにゃを見比べてみる事にします。

上が共分散の計算に指数カーネルを用いたもの、下がガウシアンカーネルを用いた物です。カーネルの選び方によって、だいぶ概形が変わるようですね。ヒートマップは分散共分散行列を可視化したものです。ヒートマップからは、隣接するデータの相関が強く、離れるほど相関が弱くなる事が分かります。この他にもカーネルは数多く存在しており、時に複数のカーネルを組み合わせたりすることもできます。カーネルによってグラフの概形や振る舞いが変わるようですね。

2-4　SOGPまとめ

ここまでで、なんとなくSOGPの雰囲気は掴めたかなと思います。重要なポイントは、以下の2つです。

ガウス過程は関数をサンプルするもので、その正体は多変量正規分布（無限次元）
分散共分散行列の各要素はカーネルなるものを使って計算、カーネルによってグラフの概形や振る舞いが変わる。

次からMOGPを述べていきます。

3．Multi-output Gaussian Process (MOGP)

さて、ここからが本題です。以下のような場面を考えてみます。

3つの色のデータが存在していますが、これらは系列単位で互いに相関しているものとしましょう。系列間の相関を捉えながらどのようにしてガウス過程を適用していきましょうか……。結論、これは圧倒的力技です。何が力技かというと、系列を順番に愚直に並べて1つの多変量正規分布にまとめてしまうという事をします。例えば、系列が $\displaystyle P$ 個あり、それぞれの系列でデータの数が $\displaystyle N$ 個あるとします。ここで、それぞれの系列は多変量正規分布に従う（ガウス過程に従う）ものとします。MOGPでは、この出力を縦に全部並べて長さ $\displaystyle NP$ の確率変数ベクトルを作成します。そうすると、長さ $\displaystyle NP$ の期待値ベクトル、 $\displaystyle NP×NP$ 個の要素を持つ分散共分散行列をパラメータとして多変量正規分布を作ることができます*5。気持ち良いくらいに力技です。

多出力のガウス過程を考える時は色々な手法が存在していますが、この記事では主にLinear Model of Coregionalization (LMC)というものを紹介し、実際にHello Worldしていきたいと思います。LMCはIntrinsic Coregionalization Model (ICM)というものとSemiparametric Latent Factor Model (SLFM)というモデルを拡張したものとなっており、この3つは全て同じ考え方の下で構築されます。以下、順を追って述べていきたいと思います。

3-1　前提

互いに相関のある系列ですが、手元のデータは横軸の値が全ての系列で共有されているとは限りません。具体的には、以下の図で示される様子です。

右側のデータはデータがサンプルされる横軸の値がすべて共通となっていますが、左側はそうはなっていません。このようなデータはそれぞれHeterotopic dataとIsotopic dataという名前で区別されています。HeterotopicとIsotopicの区別をされることで何が変わるのかと言えば、数式の記述が若干変わります。本記事ではどっちで考えても大丈夫です。加えて、それぞれの系列で同じ数だけデータが取れているとも限らないので、そこの仮定も頭に入れる必要があります。とは言え、とりあえずは簡単のために、Isotopic dataで全ての系列で同じ数のデータが取れたことを想定しておきます。

3-2　概要

・イメージ

ICM、SLFMと来て最終的にLMCを説明していくわけですが、この3つには共通点があります。それは、どのモデルも背後に潜在的なガウス過程を仮定する事です。図を見たほうが早いですね。

詳しい記号はすぐ後で説明しますので、こういう事をしてるんだなあ程度で大丈夫です。図の $\displaystyle f$ はモデル化したい過程であり、それが $\displaystyle g$ と書かれるものによって表されることを示しています。ひとまず $\displaystyle g$ が $\displaystyle f$ の上の方にいるんだなあ的な感じでOKです。

・式の展開

ICMとSLFMはなんか似ている事が分かりました。ここから少しだけ式を入れていきます。上で示した図は系列が3つで少し多い感じがするので、 $\displaystyle {\bf{f_1, f_2}}$ の2つに絞る事にします。さらに、2つの系列でそれぞれ2点が得られている事とします。つまり、 $\displaystyle {\bf{f_1}}=(f_1(x_1)　f_1(x_2))^T, {\bf{f_2}}=(f_2(x_2)　f_2(x_2))^T$ です。MOGPの冒頭でMOGPは力技であるという事を述べたので、力技をしてみましょう。

1つの多変量正規分布にまとめてみました。平均は今まで通り全て0にしておくことにします。さて、問題は分散共分散行列になる訳ですが、どのように計算されるのかを見ていく事にします。ひとまず分散共分散行列だけを取り除くと、以下のようになります。

対角成分は分散でそれ以外の成分は共分散である事が分かります。これも最初の方に述べた多変量正規分布の時と話は全く同じです。共分散は期待値を用いて、 $\displaystyle Cov(X, Y)=E[XY] - E[X][Y]$ と計算できる訳なので、1つ1つ地道に計算していけばMOGPの分散共分散行列も計算することができます。ここまではICMからLMCまですべて共通する話です。ここから先、 $\displaystyle f$ が潜在的なガウス過程 $\displaystyle g$ を用いてどのように構築されるかによって、ICMからLMCで分散共分散行列が変わってきます。

3-3　Intrinsic Coregionalization Model (ICM)

ICMは潜在的なガウス過程は1つで、そこからのサンプルによって $\displaystyle f$ を構成します。ここでは簡単のために、サンプルを1つにして次のように $\displaystyle f$ を構築します。

$\displaystyle f_1(x_1) = \omega_1^{(1)} g^{(1)}(x_1)$
$\displaystyle ...$
$\displaystyle f_2(x_2) = \omega_2^{(1)} g^{(1)}(x_2)$

ここで、 $\displaystyle g^{(1)}$ は潜在的なガウス過程 $\displaystyle g$ からのサンプルであるとします*6。この時、共分散を計算してみましょう。共分散の公式は上で示したのでそのまま代入します。例として、 $\displaystyle Cov( f_1(x_1), f_2(x_1))$ を計算してみます。

$\displaystyle Cov( f_1(x_1), f_2(x_1)) = E[f_1(x_1)f_2(x_1)] - E[f_1(x_1)]E[f_2(x_1)]$
$\displaystyle =\omega_1^{(1)}\omega_2^{(1)} ( E[ g^{(1)}(x_1)g^{(1)}(x_2) ] - E[g^{(1)}(x_1)]E[g^{(1)}(x_2)] ) = \omega_1^{(1)}\omega_2^{(1)} k(x_1, x_2)$

$\displaystyle E[f_1(x_1)]E[f_2(x_1)]$ の項は平均0を仮定しているので現れません。最後の式変形はそのまま共分散の形が出てきたので、ここをカーネルで計算させます。これを真面目にすべて計算すると次のような分散共分散行列が出てきます。

ここで、 $\displaystyle B$ は2行2列の係数の行列で、例えば $\displaystyle b_{11} = (\omega_1^{(1)})^2, b_{21} = \omega_1^{(1)}\omega_2^{(1)}, ...$ というようになっています。 $\displaystyle K$ は2行2列のカーネルを要素に持つ行列です。行列の中身がブロックに区分されているように見える訳ですが、このような物はクロネッカー積というものですっきり書く事が出来ます。最終的に導出されるものを先に見てみましょう。

最右辺の豆電球みたいな演算がクロネッカー積です。クロネッカー積の説明はとりあえず省略しておきます。適当にwikipediaとか貼っておきます。
ja.wikipedia.org

というわけで、力技で1つにまとめた分布をすっきりした形で記述することができました。

3-4　Semiparametric Latent Factor Model (SLFM)

SLFMもICMと同様に背後に潜在的なガウス過程を仮定します。しかし、ICMと違うのは潜在的なガウス過程を複数仮定することです。例えば、次のようにして $\displaystyle f$ を構築することとします。

$\displaystyle f_1(x_1) = \omega_{1, 1} g_1(x_1) + \omega_{1, 2} g_2(x_1)$
$\displaystyle ...$
$\displaystyle f_2(x_2) = \omega_{2, 1} g_1(x_2) + \omega_{2, 2} g_2(x_2)$

ここで、 $\displaystyle g_1, g_2$ はそれぞれ異なるガウス過程です。ICMで計算したのと同じように、共分散を計算してみましょう。例えば、 $\displaystyle Cov( f_1(x_1), f_2(x_1))$ を計算してみましょうか。

$\displaystyle Cov( f_1(x_1), f_2(x_1)) = E[f_1(x_1)f_2(x_1)] - E[f_1(x_1)]E[f_2(x_1)]$
$\displaystyle = \omega_{1, 1}\omega_{2, 1} Cov( g_1(x_1), g_1(x_1) ) + \omega_{1, 2}\omega_{2, 2} Cov( g_2(x_1), g_2(x_1) )$
$\displaystyle = \omega_{1, 1}\omega_{2, 1} k_1(x_1, x_1) + \omega_{1, 2}\omega_{2, 2}k_2(x_1, x_1)$

計算の過程で $\displaystyle E[f_1(x_1)]E[f_2(x_1)]$ の項は平均0を仮定しているので現れません。加えて、 $\displaystyle Cov( g_1(x_1), g_1(x_1) ), Cov( g_2(x_1), g_2(x_1) )$ はカーネルそのものです。以後、こんな形で頑張って計算すると、最終的には以下のようになります。

ICMと比べると別の分散共分散行列が出てきましたね。

3-5　Linear Model of Coregionalization (LMC)

最後にLMCです。まずは以下の図を見てみましょう。

LMCはICMとSLFMを合わせたものとなります。SLFMでまずは潜在的なガウス過程の個数を決め、各ガウス過程から複数のサンプルを取得する…といったような感じです（複数のサンプルを取得する記述は参考文献[3]で記載されていますが、それ以外で僕の見た中でLMCを説明している論文のサンプル数は1です。この辺りどうなんですかね？）。まずは、これまでの流れ通りに $\displaystyle g$ を用いて $\displaystyle f$ を構成してみましょう。ここはいきなり一般化された定義から見ていきます。

$\displaystyle f_n(x) = \sum_{q=1}^{Q} \sum_{r=1}^{r_q} \omega_{n, q}^{(r)} g_{q}^{(r)} (x)$

ここから共分散を計算していくわけですが、1から計算する必要はありません。すでにICMとSLFMで分散共分散行列の話は済ませてあるので、それを利用して求めてしまいます。LMCでは最終的に以下のような結果となります。

あれれ？ICMでいっぱいサンプルしたとすればそれはどこで反映されるの？という感じに見えますが大丈夫です。1つのガウス過程から2つ以上サンプルしたとしても $\displaystyle B⊗K$ の形で書くことができます。こんな感じでLMCはあっさりと説明を終えることができました。

4．MOGPのHello World

ここまで非常に長くなってしまいましたが、いよいよ実験してみたいと思います。とりあえずはプログラムを動かして結果を図示してみる所で終わりたいと思います。データはMOGPの冒頭に示したグラフ……ですが、ページ戻るのも大変なので、もう一度貼っておきましょうか。

MOGPはGPyを利用して動かしました。GPyを選んだ理由は特にありません。データセットの作成を含めたプログラムの全体は1番最後に貼っておきます。ひとまず核となるプログラムだけ載せる事とします。いよいよLCMにHello Worldです。まずは以下のコードを動かしてみます。

K1 = GPy.kern.Bias(1)
K2 = GPy.kern.Linear(1)
lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K1, K2])
m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
m.optimize()

初めてBiasという名前のカーネルを聞いたんですけど、こんなのあったんですね。Biasという名前のカーネルだけだと平均辺りで横一直線しか引っ張られなかったので、Linearも加えてみました。結果は以下のようになりました。

右肩上がりや右肩下がりの傾向は掴んでいます。次はシンプルにガウシアンカーネルだけにしてみましょう。

K = GPy.kern.RBF(1)
lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K])
m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
m.optimize()

良い感じのぐにゃぐにゃ具合ですね。1つ1つの系列でSOGPを掛けるよりも良い気がします。

5．終わりに

この記事ではSOGPとMOGPについて述べ、MOGPの中でもLMCに着目してHello Worldしてみました。系列全体の相関を考慮してモデル化できるのはすごく魅力的でした。しかも、全部1つにまとめる力技で動いているとは、、、多変量正規分布さまさまですね。ちなみに、この記事では以下の内容をごっそり抜かしました。

(前準備) 2つの確率変数ベクトルの同時確率密度関数と条件付確率密度関数
(前準備) ベイズの定理
(SOGP) ノイズを含んだガウス過程
(SOGP) 予測分布の導出
(SOGP) 各種推論（変分ベイズ、期待値伝播法、MCMC…）
(MOGP) MOGPの根本的な理論

などなど……

あくまでガウス過程の雰囲気をちょこっと掴もう！というHello World目的で記事を書いたのでとりあえず良いかな？という感じがします。なんとなく雰囲気が掴めた所で、より具体的な理論を学びたいと思った時には、何かしらの書籍や論文で上の内容を学ぶと良いと思います。僕も現在進行形で勉強中です！以上、Advent Calendarの23日目でした。皆様良いクリスマスをお過ごしください！それでは！

参考文献

[1] M. van der Wilk et al., "A Framework for Interdomain and Multioutput Gaussian Processes", 2020.
[2] H. Liu, J. Cai and Y-S. Ong, "Remarks on Multi-Output Gaussian Process Remarks", 2017.
[3] M. A. A´ lvarez, L. Rosasco and N. D. Lawrence, "Kernels for Vector-Valued Functions: a Review", 2012.
[4] 持橋大地，大羽成征，”ガウス過程と機械学習"．
[5] C. M. Bishop, ”パターン認識と機械学習（上下）”．

※GPyを用いて遊んでみましたが、以下を参考にプログラムを作成しました。
nbviewer.jupyter.org
github.com

プログラム

# seibisi(2歳)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import GPy

def makedata(amp, k, phi, N):
    X = np.random.rand(N)[:,None]
    Y = amp * np.sin(X/k+phi) + np.random.rand(X.size)[:,None]
    return X, Y

def plotall(m):
    fig = plt.figure(figsize=(15, 12))
    # output 1
    ax1 = fig.add_subplot(311)
    m.plot(fixed_inputs=[(1,0)], which_data_rows=slice(0, 50), ax=ax1)
    ax1.set_title("f1(x)")
    # output 2
    ax2 = fig.add_subplot(312)
    m.plot(fixed_inputs=[(1,1)], which_data_rows=slice(50, 100), ax=ax2)
    ax2.set_title("f2(x)")
    # output 3
    ax3 = fig.add_subplot(313)
    m.plot(fixed_inputs=[(1,2)], which_data_rows=slice(100, 150), ax=ax3)
    ax3.set_title("f3(x)")

def main():   
    # 1. データセットの作成
    X1, Y1 = makedata(1, 0.2, 0, 50)
    X2, Y2 = makedata(2, 0.3, 2, 50)
    X3, Y3 = makedata(0.8, 0.15, 1.5, 50)
    # 2. 当てはめ
    K = GPy.kern.RBF(1)
    lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K])
    m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
    m.optimize()
    # 3. プロット
    plotall(m)
    
if __name__ == "__main__":
    main()

*1:去年も僕はAdvent Calendarを書きました（カーネルロジスティック回帰を自前で実装して遊んだ）

*2:日本語の情報も少ない？良い日本語訳をご存知の方は教えて下さい！

*3:というよりそのもの

*4:カーネルが要素の行列はグラム行列とか名前がついていますが、しばらくは分散共分散行列という言葉のまま進めていきたいと思います。

*5:系列の1つ1つが多変量正規分布に従うならば、それを結合したものもまた多変量正規分布に従うだろうというロジックです。

*6:複数考える場合は独立なサンプルとして扱います

2020-08-18

反復重み付き最小二乗法で遊ぶ

統計

1．はじめに

積読が増えてきてしまった*1ので、最近はDobson本を読み始めたのですが、反復重み付き最小二乗法（IRLS）なるとても興味深いものを見つけました。面白そうだったからとりあえず実装してみたというのが今回の記事です。本記事は、IRLSを実装するために最低限必要な数式しか載せていないので、文章や式の展開がだいぶ飛躍しています。これに関連して、記事中で出てくる記号について先に述べておくと、 $\displaystyle {\boldsymbol{\beta}}$ は推定したいパラメータであり、 $\displaystyle {\bf{b}}$ は $\displaystyle {\boldsymbol{\beta}}$ の推定値となります。

2．指数型分布族と一般化線形モデル

指数型分布族から述べます。指数型分布族は次の式で表されます。

$\displaystyle f(y;\theta) = exp(a(y)b(\theta)+c(\theta)+d(y))$
独立な確率変数 $\displaystyle Y_1, ..., Y_N$ が従う分布を、指数型分布族に属する分布であると仮定します。さらに、以下に示す関係を仮定します。

$\displaystyle E[Y_i] = \mu_i$
$\displaystyle g(\mu_i) = {\bf{x}}_i^T {\boldsymbol{\beta}} = \eta_i$

ここで、 $\displaystyle {\bf{x}}_i^T$ は $\displaystyle i$ 番目のデータ $\displaystyle Y_i$ に対する説明変数（p×1ベクトル）であり、 $\displaystyle {\boldsymbol{\beta}}$ は推定したいパラメータです。一般化線形モデルでは、連結関数 $\displaystyle g(\mu_i)$ が非常に重要なパーツとなります。

3．パラメータの推定（IRLS）

3.1　とりあえず結論

IRLSの本題に移りますが、冒頭にも述べた通り実装に不要な数式は一気にすっ飛ばしていくので、とりあえずIRLSの反復式を示します。反復式は次のようになります。

$\displaystyle {\bf{X}}^T {\bf{WXb}}^{(m)} = {\bf{X}}^T {\bf{Wz}}$

ここで、 $\displaystyle {\bf{X}}$ は説明変数の行列でありN×p行列です。 $\displaystyle {\bf{b}}$ はパラメータの推定値ベクトルでp×1ベクトルです。 $\displaystyle {\bf{W}}$ はN×Nの対角行列であり、 $\displaystyle (i, i)$ 成分は次の式で表されます（式中の偏微分の箇所は後で説明します）。

$\displaystyle \omega_{ii} = \frac{1}{Var(Y_i)} (\frac{\partial \mu_i}{ \partial \eta_i})^2$

最後に、 $\displaystyle {\bf{z}}$ はN×1ベクトルであり、 $\displaystyle i$ 番目の要素は次の式で表されます。

$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + (y_i - \mu_i)(\frac{\partial \eta_i}{\partial \mu_i})$

IRLSを実装する時は、 $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ を逐一計算していきます。次に、この反復式を導出するに当たって必要な流れを述べていきます。

3.2　導出

最尤推定の流れから導出していきます。指数型分布族と一般化線形モデルの部分で述べた仮定に基づいて、まずは対数尤度関数を求めていきます。 $\displaystyle i$ 番目のデータに対する対数尤度関数 $\displaystyle l_i$ は次の式で表されます。

$\displaystyle l_i = y_ib(\theta_i) + c(\theta_i) + d(y_i)$

$\displaystyle a(y_i)$ が $\displaystyle y_i$ となっているのは、一般化線形モデルの場合は正準形を持つ分布が前提となっているからです。ここからパラメータの最尤推定値を求めていくわけですが、さっくり流れだけ追っていきます。パラメータの最尤推定値を得るために、対数尤度関数を微分します。とりあえず簡単のために、 $\displaystyle i$ 番目のデータに対する対数尤度関数を微分します。

$\displaystyle \frac{\partial l_i}{\partial \beta_j} = \frac{\partial l_i}{\partial \theta_i} \frac{\partial \theta_i}{\partial \mu_i} \frac{\partial \mu_i}{\partial \beta_j}$

対数尤度関数は合成関数になっている*2ため、微分の連鎖法則を用いて微分していきます。さらに、右辺の3つ目の偏微分は次のように微分の連鎖法則で微分できます。

$\displaystyle \frac{\partial \mu_i}{\partial \beta_j} = \frac{\partial \mu_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta_j} = \frac{\partial \mu_i}{\partial \eta_i}x_{ij}$

一番最後の式の偏微分はとりあえず結論で述べた所にも現れています。この偏微分は（記事には載せませんが）スコアや情報行列にも現れてきており、解析的に反復の式を導出する際はこれを計算する必要があります。

あとは偏微分から良い感じになんやかんやしてなんやかんやする*3と、IRLSの反復計算式を求める事が出来ます。途中経過は省略です。

3.3　例．ロジスティック回帰

実験ではロジスティック回帰で遊んでいきます。Dobson本では第7章で触れられています*4が、スコアと情報行列を出して終わっているので、IRLSで計算すべき $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ までは述べてくれていません。従って、自力で導出していきます。計算が間違えている可能性もあるので、あらかじめごめんなさいしておきます。

ひとまず重要なのは、次の2つです。

$\displaystyle \omega_{ii} = \frac{1}{Var(Y_i)} (\frac{\partial \mu_i}{ \partial \eta_i})^2$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + (y_i - \mu_i)(\frac{\partial \eta_i}{\partial \mu_i})$

ここで、 $\displaystyle Y_i～binomial(n_i, \pi_i)$ であり、ロジスティック回帰では期待値と連結関数は次のようになります。

$\displaystyle E[Y_i] = \mu_i = n_i \pi_i$
$\displaystyle g(\pi_i) = \eta_i = ln \frac{\pi_i}{1-\pi_i} = {\bf{x}}_i^T {\boldsymbol{\beta}}$

$\displaystyle w_{ii}$ 中の $\displaystyle Var(Y_i)$ は $\displaystyle n_i \pi_i (1-\pi_i)$ なのでとりあえずOKです。次に、どちらにも共通して現れる偏微分を求めていきます。偏微分の結果は次のようになります。

$\displaystyle \frac{\partial \mu_i}{\partial \eta_i} = \frac{\partial (n_i \pi_i) }{\partial \eta_i} = n_i \frac{\partial \pi_i}{\partial \eta_i} = n_i (\frac{\partial \eta_i}{\partial \pi_i})^{-1} = n_i \pi_i (1-\pi_i)$

第3式から第4式は逆関数の微分を利用しました。これで、 $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ は次のようにすっきり求める事が出来ました。

$\displaystyle \omega_{ii} = n_i \pi_i (1-\pi_i)$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + \frac{y_i - n_i \pi_i }{ n_i \pi_i (1-\pi_i) }$

これで完成で、あとは反復式に入れるだけです。

4．実験

4.1　IRLSの実装

さて、長々と述べてきたわけですが、IRLSの実装に重要なポイントをまとめていきます。

反復式

$\displaystyle {\bf{X}}^T {\bf{WXb}}^{(m)} = {\bf{X}}^T {\bf{Wz}}$
反復式では、 $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ を計算する必要があります。

ロジスティック回帰のIRLS

ロジスティック回帰では、 $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ は次の式で表す事が出来ました。

$\displaystyle \omega_{ii} = n_i \pi_i (1-\pi_i)$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + \frac{y_i - n_i \pi_i }{ n_i \pi_i (1-\pi_i) }$
あとは上手く反復の手順を考えるだけです。

反復の手順

パラメータが多いと可視化に困るので、 $\displaystyle {\bf{x}}_i$ と $\displaystyle {\bf{b}}$ を次のように決めておきます。

$\displaystyle {\bf{x}}_i^T = [ 1　x_i], {\bf{b}}^T = [ b_1　b_2]$

$\displaystyle m$ ステップ目の推定したパラメータを $\displaystyle \hat{\bf{b}}^{(m)}$ とした時、反復の手順を次に示します。

手順1：初期値を決める→パラメータの初期値を $\displaystyle \hat{ \bf{b} }^{(0)}$ 、確率の最尤推定値の初期値を $\displaystyle {\hat{\pi_i}}^{(0)}$ として与えます。
手順2： $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ を計算
手順3：反復式を次のように変形して、パラメータを推定

$\displaystyle \hat{ \bf{b} }^{(1)} = ({\bf{X}}^T {\bf{WX}} )^{-1} {\bf{X}}^T {\bf{Wz}}$
手順4：パラメータの推定後、 $\displaystyle {\hat{\pi_i}}^{(1)}$ を次の式で計算
$\displaystyle {\hat{\pi_i}}^{(1)} = \frac{exp({\hat{b_1}}^{(1)} + \hat{b_2}^{(1)} x_i )}{1 + exp({\hat{b_1}}^{(1)} + \hat{b_2}^{(1)} x_i ) }$
手順5：手順1の $\displaystyle \hat{ \bf{b} }^{(0)}, {\hat{\pi_i}}^{(0)}$ を $\displaystyle \hat{ \bf{b} }^{(1)}, {\hat{\pi_i}}^{(1)}$ に変更して、手順2から繰り返す

これをパラメータが収束するまで繰り返していきます。

4.2　実験

まず初めにデータセットを作成します。 $\displaystyle \beta_1 = 2.0, \beta_2 = 3.0$ とし、 $\displaystyle (x_i, \pi_i)$ を20個作成しました。この時、生成したデータセットに対してロジスティック回帰モデルを当てはめ、IRLSでパラメータを求める様子をアニメーション化してみました。

左はフィットの様子で、紫のデータは作成したデータセットです。右はパラメータの推移の様子で、青があらかじめ設定したパラメータです。左右共通して赤が推定した結果です。IRLSはあっという間に収束してしまうようで、じわじわパラメータが推定される様子が見えたら良かったのですが、なんだかつまらない感じになってしまいました*5。最終的なパラメータの推定値は、 $\displaystyle b_1 = 1.99, b_2 = 3.01$ でした。

5．まとめ

ちゃんとうごかすことができてとてもたのしかったです

プログラム

こっそりあとがきですが、実はこの記事を書いた段階ではDobson本は第4章までしか読めておらず、収束の判定まではプログラムに実装できませんでした*6。安直にfor文で回しただけのプログラムになってしまいましたが、IRLSのコアとなる部分はfor文の中の4行です。もう少し上手い事プログラムを組めばより短くなるかもしれません。

# seibisi
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.random import *
from matplotlib import animation

def make_p(beta1, beta2, x):
    return np.exp(beta1 + beta2 * x) / (1 + np.exp(beta1 + beta2 * x))
    
def main():
    # パラメータ
    beta = [2.0, 3.0]
    x = np.linspace(-5, 5, 20)
    pi = make_p(beta[0], beta[1], x)
    # データセットの作成
    n = [round(normal(50, 5)) for i in range(len(x))]
    y = [round(nn * pp) for nn, pp in zip(n, pi)]
    # 初期値
    beta_hat = [0, 0]
    pi_hat = [0.5 for i in range(len(y))]
    X = np.array([np.ones(len(y)), x]).T
    # IRLS
    for i in range(10):
        # 1. zとWの計算
        z = [beta_hat[0] + beta_hat[1] * xx + (yy - nn * pp)/(nn * pp * (1-pp)) for xx, yy, nn, pp in zip(x, y, n, pi_hat)]
        W = np.diag([nn * pp * (1-pp) for nn, pp in zip(n, pi_hat)])
        # 2. パラメータの計算
        beta_hat = np.dot(np.linalg.inv(np.dot(np.dot(X.T, W), X)), np.dot(np.dot(X.T, W), z))
        # 3. 予測
        pi_hat = make_p(beta_hat[0], beta_hat[1], x)

*1:博論を書いているうちに積んでしまっていた

*2:えぐい

*3:ここからスコアと情報行列を導き、情報行列を行列の積で表すと出てきます

*4:実は第4章で読むのが止まっている

*5:7, 8ステップ目で完全に収束しました。pythonかRかなんかを使って収束するまでのステップ数を比較すればもっと面白かったかもしれない

*6:最後まで読んでから実装しよう

2020-06-11

標本自己相関関数で遊ぶ

統計

はじめに

（またまた）久しぶりの更新となりました。最近は研究で時系列を扱う機会があり、時系列の勉強をしています。とりあえず北川源四郎「時系列解析入門」（以下、北川時系列）と沖本竜義「経済・ファイナンスデータの計量時系列分析」（以下、沖本時系列）の2冊を手に入れ、まずは北川時系列を読み始めました。とりあえず勉強するがてら手も動かそうという事で、標本自己相関関数で遊んでみたいと思います。

定義

定常な時系列 $\displaystyle {\bf y}=\{ y_1, y_2, ..., y_T\}$ が与えられた時、標本自己共分散関数 $\displaystyle {\hat C_k}$ と標本自己相関関数 $\displaystyle {\hat R_k}$ は次の式で定義されます。

$\displaystyle {\hat C_k} = \frac{1}{N} \sum_{t=k+1}^{T} (y_t - {\hat \mu}) (y_{t-k} - {\hat \mu} )$
$\displaystyle {\hat R_k} = \frac{ {\hat C_k} }{ {\hat C_0} }$
定義を見ればああなるほどという感じですが、「標本」自己共分散『関数』という事で関数自体が標本？？になったので、実験して遊んでみたいと思います。きっと時系列のぐにゃぐにゃ1本1本が標本であり、それに対する自己相関関数だから「標本自己相関関数」なのかな？と思います（合ってるのか自身が無い）。

実験

ホワイトノイズとランダムウォークをそれぞれ100本作成し、それぞれで標本自己相関関数をプロットしてみます。

図の上段がホワイトノイズ、下段がランダムウォークです。左側が時系列であり、右側が自己相関関数です。ホワイトノイズはどのぐにゃぐにゃであっても標本自己相関関数は割と同じような感じになっていますが、ランダムウォークはだいぶ異なっているみたいです。全然話は変わりますが、ランダムウォークのグラフがNintendo 64のポケモンスタジアムでポケモンが破壊光線を打った時のグラフィックっぽくて懐かしさを覚えました。初代の破壊光線は凶悪でしたね。

終わりに

とても短い記事となってしまいましたが、今回は標本自己相関関数で遊んでみました。独立性が満たされないデータの分析は非常に興味があったので、とても面白く勉強できています。少しずつ読み進めて、面白いのがあったらまた手を動かして遊んでみたいと思います。

プログラム

# seibisi
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def main():
    N = 500
    lag = 200
    y1 = []
    y2 = []
    autocorr1 = []
    autocorr2 = []
    for i in range(100):
        # ホワイトノイズとランダムウォーク
        yy1 = pd.Series(np.random.normal(0, 1, N)) 
        yy2 = np.cumsum(yy1)
        # 自己相関
        autocorr1.append([yy1.autocorr(l) for l in range(lag)])
        autocorr2.append([yy2.autocorr(l) for l in range(lag)])
        y1.append(yy1)
        y2.append(yy2)
    # 描画
    t = [i for i in range(N)]
    k = [i for i in range(lag)]
    """
    for y in y1:
        plt.plot(t, y)
    plt.xlabel("t")
    plt.ylabel("y")
    plt.show()
    """
    """
    for ac in autocorr2:
        plt.plot(k, ac)
    plt.xlabel("lag")
    plt.ylabel("corr") 
    plt.show()
    """
    
if __name__ == "__main__":
    main()

2020-05-04

Maternカーネルで遊ぶ

機械学習

はじめに

久々の更新となりました。最近、研究でガウス過程回帰を用いる用があり、MLPシリーズのガウス過程と機械学習で勉強していた所、Maternカーネルなる面白い物を見つけました。Maternカーネルはとあるパラメータ（後述）を変える事で、指数カーネルになったりRBFカーネルになったりするようです。非常に面白いですよね。今回はこのMaternカーネルなる不思議なカーネルで遊んでいきたいと思います*1 *2。カーネルとはなんぞ？？みたいな話はしないので、ある程度カーネルの話について知っていると読みやすいと思います。

Maternカーネル

Maternカーネルの定義

まず初めに式を示します。

$\displaystyle k_{\nu} ( {\bf x, x^{'} } ) = \frac{ 2^{1-\nu} }{ \Gamma (\nu) } ( \frac{ \sqrt{2\nu}r }{\theta} )^{\nu} K_{\nu} ( \frac{ \sqrt{2\nu}r }{\theta} )$
ここで、 $\displaystyle K_{\nu}$ は第2種の変形ベッセル関数、 $\displaystyle \theta$ はスケールパラメータ、 $\displaystyle r = | {\bf x - x^{'} } |$ です。

色々なMaternカーネル

定義に現れる $\displaystyle \nu$ が関数の滑らかさと関わりがあり、 $\displaystyle \nu = \frac{1}{2}$ の時は指数カーネル、 $\displaystyle \nu = \infty$ の時はRBFカーネルと一致します。他にも、Matern3, Matern5なども存在しているようです。ガウス過程と機械学習では、Matern3とMatern5の図が現れており、ああ確かにパラメータ変えれば滑らかになるんだなあという気持ちになります。

他のパラメータは使っちゃダメなの？

というのが今回の遊ぶポイントです。ある時は指数カーネル、またある時はRBFカーネルなど、パラメータによって滑らかさが変化してくれるのは良いけど、どうしてそのパラメータじゃなきゃダメなのでしょうか*3。また、パラメータを規則的に変化させると、滑らかさも規則的に変化するのでしょうか？

ちなみに、ガウス過程はGPyやscikit-learnなど、有名どころのパッケージを入れる事で簡単に遊ぶことができます。例えば、scikit-learnでMaternカーネルを調べると、次のような説明が載っています。

The parameter nu controlling the smoothness of the learned function.
（省略）
Note that values of nu not in [0.5, 1.5, 2.5, inf] incur a considerably higher computational cost (appr. 10 times higher) since they require to evaluate the modified Bessel function. Furthermore, in contrast to l, nu is kept fixed to its initial value and not optimized.

scikit-learn.org

指数カーネル、Mattern3, Mattern5、RBFカーネル以外を指定すると、やたらと計算量食うから気をつけろよとの記述があります。でもそんなことは関係ありません。早速遊んでみましょう。

実験

仮想データセットの作成

まずはデータセットの作成です。次の式からデータセットを作成します。

$y = |sin2x| + \epsilon, \epsilon～N(0, 0.1^2)$

実験では、上の式から30個のサンプルを取得し、このデータに対してガウス過程回帰を掛け、期待値と1σ範囲をプロットします。また、滑らかさを調節するパラメータ $\displaystyle \nu$ は0.01から1.5まで0.01刻みで変化させることにします。

結果

結果はこんな感じになりました。赤が予測、黄色が1σ範囲です。

最初は全然データに対応できていない上に結構変化が激しいけど、時間が経つとなんだか変化がスローモーションっぽくなりますね。規則的にパラメータを変化させても滑らかさの変化は規則的ではないような気がしますね。加えて、Mattern0 ( $\displaystyle \nu=0$ )みたいなカーネルがあったとすれば、きっと使い物にならないんだろうなあという事も副次的に分かりました。

終わりに

今回はMaternカーネルで遊んでみました。正直記事として耐える内容なのかそうでないのかだいぶ迷いましたが、GWなので（？）とりあえず書きました。ところでMaternのeの上に「ちょん」*4がついているんですけど、どうやってキーボード入力するんですかね。いちいちコピペするのも面倒だし…まあいいや

プログラム

一応載せておきます

# seibisi

from sklearn.gaussian_process import kernels
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from matplotlib import animation

def make(x):
    return np.abs(np.sin(2*x)) + np.random.normal(0, 0.1, len(x))

# 当てはめ，予測
def fit_pred(x, nu, df):
    # 当てはめ
    kern = kernels.Matern(nu=nu)
    clf = GaussianProcessRegressor(kernel=kern, alpha=1e-10, n_restarts_optimizer=30, normalize_y=True)
    clf.fit(df["x"].values.reshape(-1, 1), df["y"])
    # 予測
    exp, std = clf.predict(x.reshape(-1, 1), return_std = True)
    return exp, std
    
def main():
    # 1. データセットの作成
    x = np.linspace(-2, 2, 30)
    df = pd.DataFrame(np.transpose([x, make(x)]), columns=["x","y"])
    df = df.sort_values("x").reset_index(drop=True)
    # 2. 当てはめと予測
    nus = np.linspace(0.01, 1.5, 150)
    x_test = np.linspace(-3, 3, 200)
    exp = [None for i in range(len(nus))]
    std = [None for i in range(len(nus))]
    for i in range(len(nus)):
        exp[i], std[i] = fit_pred(x_test, nus[i], df)
    # 3. アニメーションで出力
    k = 1
    fig = plt.figure()
    axes = fig.add_subplot(111)
    ims = []
    for i in range(len(exp)):
        img = [axes.scatter(df["x"], df["y"], marker="x", color="purple")]
        img += axes.plot(x_test, exp[i], color="red")
        img += axes.plot(x_test, exp[i]-k*std[i], color="orange")
        img += axes.plot(x_test, exp[i]+k*std[i], color="orange")
        #img += axes.fill_between(x_test, exp[i]-k*std[i], exp[i]+k*std[i], color="orange", alpha=0.1)
        ims.append(img)
    ani = animation.ArtistAnimation(fig, ims, interval = 80)
    ani.save("ani.gif", writer="pillow")
    
if __name__ == "__main__":
    main()

*1:実は既に遊んでいてTwitterにも載せていたけど、敢えて記事にするほどの内容かどうかで迷っていてなんやかんやで時間が過ぎてしまった

*2:とりあえず書くことにした

*3:指数カーネルとRBFカーネルに一致する話については、ダメというよりも一致するから便利よ位の話だと思うので別に良いけど、Matern3とかMater5とか…

*4:ドイツ語のウムラウト的な…？

2020-03-15

WAICで遊ぶ

統計

1．はじめに

WAICはモデル選択のための評価基準の1つです。WAICは汎化損失を推定する量として非常に有用性の高い評価基準となっています。この辺りの理論や細かい話については、言わずと知れたベイズ統計の名著であるベイズ統計の理論と方法（渡辺ベイズ）や種々の記事に譲ること*1にして、本記事では実際に手を動かしてWAICで遊んでみたいと思います。ある程度渡辺ベイズを読まれていると、何をして遊んでいるのかが分かるかもしれません。

※第3章までは泣く泣く読んでいきましたが、第4章になってから途端に読めなくなってしまったので、とりあえず一旦整理のために手を動かしてみました…間違いや解釈違いは多々あると思います…

2．汎化損失とWAIC

さて、実際に遊んでいく前に、少しだけ汎化損失とWAIC周りの話を簡単に述べていきます。初めに、汎化損失 $\displaystyle G_n$ は以下の式で定義されます。

$\displaystyle G_n = - \int q(x)logp^*(x)dx$

汎化損失は真の分布 $\displaystyle q(x)$ に対して、予測分布 $\displaystyle p^*(x)$ がどの程度推測できているかを量的に示しています。

しかし、実際には真の分布 $\displaystyle q(x)$ は未知であるため、正確な汎化損失を計算することはできません。従って、何らかの形で「推定」が必要になります。このような状況で、WAICは汎化損失の推定量として大きな役割を果たします。WAICは次の式で定義されます。

$\displaystyle W_n = T_n + \frac{\beta V_n}{n}$

ここで、 $\displaystyle T_n$ は経験損失、 $\displaystyle V_n$ は汎関数分散、 $\displaystyle \beta$ は逆温度、 $\displaystyle n$ はデータの数です。経験損失と汎関数分散については以下の式で計算されます。

$\displaystyle T_n = - \frac{1}{n} \sum_{i=1}^{n} log E_{\omega} [ p(X_i | \omega) ]$
$\displaystyle V_n = \sum_{i=1}^{n} \{ E_{\omega} [ ( logp(X_i | \omega) )^2 ] - E_{\omega} [ logp(X_i | \omega) ]^2 \}$

ここで、 $\displaystyle \omega$ は推定されたパラメータ、 $\displaystyle E_{\omega} []$ は推定された事後分布で期待値を取る操作になります。経験損失も汎関数分散も、パラメータさえ推定してしまえば計算できる量なので、WAICは手元にあるデータから完全に計算することができます。

加えて、汎化損失とWAICについて、次に示す非常に興味深い関係が存在しています。

$\displaystyle E[G_n] = E[W_n] + o(\frac{1}{n})$

つまり、データを取ってはWAICを計算して、またデータを取ってはWAICを計算して…を繰り返して期待値を取ると、汎化損失の期待値に一致するというものです。汎化損失は未知の分布 $\displaystyle q(x)$ が入っていることで実際に計算することができない量でしたが、平均的に考えればこの量を見事に推定できているという関係です。非常に面白い関係です。

WAICを利用してモデル選択を行う際は、WAICが小さいモデルを選択すると良いわけですが、WAICは確率変数であることを頭に入れておくことが大切です。

3．シミュレーション

というわけで、早速WAICで遊んでみましょう。

未知の分布 $\displaystyle q(x)$ の作成

今回のシミュレーションでは、未知の分布を混合正規分布としました。具体的には、混合数を3（混合比＝ $\displaystyle (0.2, 0.5, 0.3)$ ）とし、 $\displaystyle N(0, 1), N(2, 1), N(4, 1)$ を混ぜています。この $\displaystyle q(x)$ に対して、①通常の正規分布②混合数2の混合正規分布③混合数3の混合正規分布の3パターンを確率モデルとして仮定したときに、それぞれのWAICがどのようにふるまうのかを見ていきます。真の分布は次の図で示す通りです。

実験結果1　WAICのHello World

サンプリングするデータ数を50にした時のヒストグラムをはじめに示します。左側の赤いグラフが $\displaystyle q(x)$ 、右側が実際のデータです。

このデータに基づいて、予測分布 $\displaystyle p^*(x)$ を求めていきます。今回はシミュレーションなので、すべての分散を既知（分散=1ですべての正規分布に渡って等分散）し、すべての平均の事前分布を $\displaystyle N(0, 1)$ としておきました。

このデータに基づいて、作成された予測分布を次に示します。左から、通常の正規分布、混合数2、混合数3の予測分布です。

汎化損失はそれぞれ2.686、2.000、1.984となり、WAICはそれぞれ2.420、2.110、2.069となりました。汎化損失とWAICはだいたい近い値を取り、モデルは混合数3の混合正規分布を選ぶと良さそうです*2

実験結果2　汎化損失とWAICの分布

汎化損失とWAICの期待値を実際に比較してみましょう。サンプリングするデータ数と事前分布を実験結果1と同じにしておき、混合数2の混合正規分布を確率モデルとして用いた*3時の汎化損失とWAICのヒストグラムを以下に示します。ここで、サンプリングからWAIC、汎化損失の計算は100回行いました。

f:id:seibibibi:20200315164357p:plain

それぞれの平均は汎化損失が2.004、WAICが2.044でした。ほとんど一致しています。すごい。

実験結果3　サンプルサイズを変更したときのWAICの分布

計算回数、確率モデル、事前分布は実験結果2と同じに揃え、データ数を50, 100, 150と変えた時のWAICの分布を示します。

f:id:seibibibi:20200315113145p:plain

データ数が増えるほど、ある1点に集中する傾向があります。面白いなあ。

4．まとめ

本記事ではWAICを使って遊んでみました。やはり数式を追うだけではなく実際に手を動かしてどうなるのかを確かめるのは楽しいですね。非常に勉強になりました。正則理論までしか読めていませんが、数学をちょこちょこ勉強して、一般理論以降の章を読めるように頑張っていきたいところです。

プログラム

WAICの実装自体はRで行い、次の記事を参考にさせていただきました。
statmodeling.hatenablog.com

パラメータのサンプリング自体はStanを利用しました。Stanのプログラムは以下になります。基本的にはStan User's Guideを参考にしましたが、その通りの実装してみると案の定上手く行かなかった*4（あるある）ので、分散を既知にした状態でプログラムを書きました。Stan User's Guideは以下になります。

https://mc-stan.org/docs/2_19/stan-users-guide-2_19.pdf

data{
  int<lower=1> K;
  int<lower=1> N;
  real y[N];
}
parameters{
  simplex[K] pi;
  ordered[K] mu;
}
model{
  vector[K] log_pi = log(pi);
  mu ~ normal(0, 1);
  for(n in 1:N){
    vector[K] lps = log_pi;
    for(k in 1:K)
      lps[k] += normal_lpdf(y[n] | mu[k], 1);
    target += log_sum_exp(lps);
  }
}
generated quantities{
  vector[N] log_likelihood;
  int index;
  real y_pred;
  for(n in 1:N){
    vector[K] lp;
    for(k in 1:K)
      lp[k] = log(pi[k]) + normal_lpdf(y[n] | mu[k], 1);
    log_likelihood[n] = log_sum_exp(lp);
  }
  index = categorical_rng(pi);
  y_pred = normal_rng(mu[categorical_rng(pi)], 1);
}

*1:譲るとは言ったものの、渡辺ベイズの中で展開される数式を追い切れていないので、完全に丸投げの形になってしまいますが…

*2:ちなみに、WAICの差がどの程度本質的なのかは僕はよくわかっていません

*3:混合数3を採用すると遅いので混合数2にしました。ただせっかちなだけです…

*4:具体的には、収束に関する警告がどっさり出てきます。ここの問題はまだ未解決なので色々調べてみたいと思います。

0．はじめに

1．前準備

2．Single-output Gaussian Process (SOGP)

2-1 SOGPの雰囲気

2-2 SOGPの定義

2-3 SOGPの分散共分散行列

2-4 SOGPまとめ

3．Multi-output Gaussian Process (MOGP)

3-1 前提

3-2 概要

・イメージ

・式の展開

3-3 Intrinsic Coregionalization Model (ICM)

3-4 Semiparametric Latent Factor Model (SLFM)

3-5 Linear Model of Coregionalization (LMC)

4．MOGPのHello World

5．終わりに

参考文献

プログラム

1．はじめに

2．指数型分布族と一般化線形モデル

3．パラメータの推定（IRLS）

3.1 とりあえず結論

3.2 導出

3.3 例．ロジスティック回帰

4．実験

4.1 IRLSの実装

反復式

ロジスティック回帰のIRLS

反復の手順

4.2 実験

5．まとめ

プログラム

はじめに

定義

実験

終わりに

プログラム

はじめに

Maternカーネル

Maternカーネルの定義

色々なMaternカーネル

他のパラメータは使っちゃダメなの？

実験

仮想データセットの作成

結果

終わりに

プログラム

1．はじめに

2．汎化損失とWAIC

3．シミュレーション

未知の分布の作成

実験結果1 WAICのHello World

実験結果2 汎化損失とWAICの分布

実験結果3 サンプルサイズを変更したときのWAICの分布

4．まとめ

プログラム

2-1　SOGPの雰囲気

2-2　SOGPの定義

2-3　SOGPの分散共分散行列

2-4　SOGPまとめ

3-1　前提

3-2　概要

3-3　Intrinsic Coregionalization Model (ICM)

3-4　Semiparametric Latent Factor Model (SLFM)

3-5　Linear Model of Coregionalization (LMC)

3.1　とりあえず結論

3.2　導出

3.3　例．ロジスティック回帰

4.1　IRLSの実装

4.2　実験

未知の分布 $\displaystyle q(x)$ の作成

実験結果1　WAICのHello World

実験結果2　汎化損失とWAICの分布

実験結果3　サンプルサイズを変更したときのWAICの分布