HELLO CYBERNETICS

深層学習、機械学習、強化学習、信号処理、制御工学、量子計算などをテーマに扱っていきます

確率モデリングと事後分布、事前分布、超パラメータ

 

 

follow us in feedly

はじめに

下記の記事に続き、あるいは知識を前提として進めていきます。

www.hellocybernetics.tech

上記の記事での内容の概要は

  • 同時分布の設計(モデリング)
  • 観測変数による条件付き分布の導出(乗法定理を利用)
  • グラフィカルモデルによるモデリング

です。これらを既に知っている人は上記の記事を参考程度に眺めるに留め、今回の記事を読んでも大丈夫だと思われます。

例題:正規分布に従っているように見えるデータ $X$ の推論

モデリング再考

さて、今回は手始めに正規分布を使った推論を確率モデリングの流れに沿って一から実施してみましょう。 まず、正規分布に従っているであろう確率変数 $X$ を考えれば、それと同時に正規分布の形状を決定する変数として$ \mu, \sigma$ を列挙することができます。ここで、$X$ が観測変数、すなわちデータに対応し、$\mu, \sigma$ が観測できない確率変数(データから推論されるもの)になります。 多くのケースではこれらをパラメータと呼ぶのですが、今は気にせず、確率モデリングの流れに沿って進めてみましょう。

同時分布に関しては

$$ p(X, \mu, \sigma) = p(\mu)p(\sigma)p(X \mid \mu, \sigma) $$

と仮定してみることにします。これは下記のようなグラフィカルモデルを考えていることに相当します。

f:id:s0sem0y:20200814184021p:plain

ここで、$X$ は正規分布に従っているように見えるのですから $p(X \mid \mu, \sigma) = \mathcal N(X \mid \mu, \sigma)$ であると考えていることになります(これは後ほど反映します)。 さて、観測変数 $X$ から $\mu, \sigma$ の条件付き分布を求めようと試みるならば乗法定理(ベイズの定理)を用いて

$$ p(\mu, \sigma \mid X) = \frac{p(X, \mu, \sigma)}{p(X)} = \frac{p(\mu)p(\sigma)p(X\mid \mu, \sigma)}{p(X)} $$

と書き下すことになります。こうして観測データに条件付けられた観測できない変数の条件付き分布が得られたことになります。

事後分布を書き下す

さて、これは数式を乗法定理と仮定した同時分布によって弄くって整理したに過ぎません。 

$$ p(\mu, \sigma \mid X) = \frac{p(\mu)p(\sigma) p(X\mid \mu, \sigma)}{p(X)} $$

ところで、$X$ が正規分布に従っているように見えるのは、$X$ を複数個、例えばここでは $N$ 個観測してからのはずです。このことを反映するために、一旦、 $X$ に関して実は $X _ 1, X _ 2, X _ 3, ..., X _ N$ と複数個の別々の確率変数が観測されたものであると考えてみることにします。すると上記の式は

$$ p(\mu, \sigma \mid X _ 1, X _ 2, X _ 3, ..., X _ N) = \frac{p(\mu)p(\sigma) p(X _ 1, X _ 2, X _ 3, ..., X _ N\mid \mu, \sigma)}{p(X _ 1, X _ 2, X _ 3, ..., X _ N)} $$

と書き換えることができます。難しいことはありません。 $X$ を $X _ 1, X _ 2, X _ 3, ..., X _ N$ に置き換えただけです。これら複数個の確率変数$X _ 1, X _ 2, X _ 3, ..., X _ N$を観測したことによって、 $\mu, \sigma$ の条件付き分布を書き下したということになります。

そして、再び個々の $X _ 1, X _ 2, X _ 3, ..., X _ N$ に着目し、それら1つ1つの確率変数 $ X _ i$ が同じ分布から独立に(他の $X _ j \ \ (j \neq i)$ の結果に無関係に)同じ分布から値が得られたのであるとすれば、独立性から

$$ p(X _ 1, X _ 2, X _ 3, ..., X _ N \mid \mu, \sigma) = p(X _ 1 \mid \mu, \sigma)...p(X _ N \mid \mu, \sigma) = \prod _ {i = 1} ^ {N} p(X _ i\mid \mu, \sigma) $$

と書き下すことができます。グラフィカルモデルではこのことを下記のように記します。

f:id:s0sem0y:20200814185637p:plain

単に確率変数 $X$ 1つを考えていたものを $N$ 個観測しているのだから、グラフィカルモデルを$N$個書いて、独立同分布($ \mu, \sigma$ は共通)であることと $X$ を繰り返し書いている部分は略記できるので、まとめただけです。 以上を踏まえ、観測データ $X ^ N = X _ 1, X _ 2, X _ 3, ..., X _ N$ (文字の方も略記 $X ^ N$ を使わせてもらいます!)による、$\mu, \sigma$ の条件付き分布は

$$ p(\mu, \sigma \mid X ^ N) = \frac{p(\mu)p(\sigma)\prod _ {i = 1} ^ N p(X _ i\mid \mu, \sigma)}{\prod _ {i = 1} ^ p(X ^ N)} $$

と書けることになります。ちなみに、ここで分子に現れている観測データ $X _ i$ に関する確率を表現している部分 $\prod _ {i = 1} ^ N p(X _ i\mid \mu, \sigma)$ を尤度関数と呼びます。 $X$ は正規分布に従っているように見えているということでしたので、このことを尤度関数に反映すれば

$$ p(\mu, \sigma \mid X ^ N) = \frac{p(\mu)p(\sigma)\prod _ {i = 1} ^ N \mathcal N(X _ i\mid \mu, \sigma)}{ p(X ^ N)} $$

と書き下すことができます。このように観測データ $X ^ N$ による条件付き分布のことを事後分布と呼びます。 さて、ここで微妙に分子にそのまま残っている $p(\mu)$ と $p(\sigma)$ を事前分布と呼びます。事前分布はグラフィカルモデルで有向グラフを書いた時に矢印が入ってきていない(出発点にしかなっていない)確率変数全てに考慮する必要があります。ここまではハッキリと述べてきませんでしたが、同時分布の設計(モデリング)とは、確率変数間のつながり(グラフィカルモデルでどこにどのように線をつなぐのか)を考えることに加えて、観測変数がどのような分布に見えているのか(今回のケースは正規分布)を尤度関数として仮定し、更に事前分布にも具体的な分布を与えることによって初めて達成されます。

事前分布と超パラメータ

さて、事前分布の設定には "観測データの分布の仮定と同様に任意性"があります。ここでは全体像と表記の説明を簡単に進めるため、事前分布として一様分布 $p(\mu) = \mathcal U(\mu \mid a, b)$ と $p(\sigma) = \mathcal U(\sigma \mid c, d)$ を与えることにします。ここで、$\mathcal U(a, b)$ は区間 $[a, b)$ に定数 $1 / (a-b)$ の値を取る確率分布です。今回 $\sigma$ は正の値であるはずなので $0<c<d$ となっていることとします。

ここで現れた $a, b, c, d$ などの事前分布の形状を決めている決め打ちの設計パラメータを超パラメータと呼びます。

ところで、こんないい加減に事前分布を決めて良いのだろうか。当然、しっかりデータのドメイン知識を反映すべきであるし、現実の複雑な問題ではドメイン知識を凝らしても良い事前分布を設定するのは簡単ではない。従って、事前分布にどのようなものを選択すべきであるかを数値的に評価する方法が必要になるし、そのような方法は用意されている。例えば今回は一様分布を選択したが、一様分布にしてもその値を取る区間の設定をどうするかという話がある。つまりハイパーパラメータの選択をどうするかという話である。これは経験ベイズ法などによって定量的に評価できる。今回はここには立ち入らない。そもそも、今回の例題のような正規分布のケースでは、事前分布の設計がどうのこうのなんてことが問題にならない場合がほとんどだろうし、多くのケースでこれほどシンプルならば、最尤推定してしまうのが手っ取り早い。

さて、超パラメータを設定しているということを明示するために

$$ p(\mu, \sigma \mid X ^ N, a, b, c, d) = \frac{\mathcal U(\mu\mid a, b)\mathcal U(\sigma\mid c, d)\prod _ {i = 1} ^ N \mathcal N(X _ i\mid \mu, \sigma)}{ p(X ^ N)} $$

などと書かれることがあります(そもそも決め打ちの定数なら、明記してもその後使われなかったりするのだが、ハイパーパラメータを後に評価する場合には文献の中でこういう書かれ方がすることがある)。また、このことをグラフィカルモデルに反映すると下記のようになります。

f:id:s0sem0y:20200814194315p:plain

グラフィカルモデルの標準的な表記

さて、グラフィカルモデルをモデリングで柔軟に使っていくためには数式とコトバによる注釈をいつでもセットにしなければならないのは面倒です。 今回の場合上記の図を見て、説明を聞いてきた人ならば $X$ が観測変数であるということを察することはできるでしょう。しかし数式よりも先にグラフィカルモデルでどういうことを仮定しているか(何が観測で、どのような関係性があるのか)をパッと表現できたら便利なはずです。そこで、標準的には観測変数と観測できない変数(通常は潜在変数と呼ぶ)をグラフィカルモデルで明確に分けるためには、観測変数はグレーで塗りつぶされ下記のような表記になります。

f:id:s0sem0y:20200814200647p:plain

例題:時系列モデル

さて、ここまでを一通り頭に入れておくと、数式とグラフィカルモデルの対応がある程度は取れるようになっており、モデリングの手助けに小さな地図が手に入ったのではないかと思います。一度、グラフィカルモデルと数式の対応が取れるようになると、数式だけではやけに複雑に見えたモデルも、頭の中で考えやすくなっているのではないでしょうか。そのことを確認すべく、時系列モデルについて少し触れていきたいと思います。

時系列データのグラフィカルモデル

さて、まずは時系列データのグラフィカルモデルを見てみましょう。このときに理解の助けになるのは、先程まで見てきた正規分布の例に戻り、観測変数 $X$ が単独で下記のグラフィカルモデルで得られているケースから丁寧に見ていくことです。

f:id:s0sem0y:20200814184021p:plain

このモデルにおいて、観測データが $N$ 個得られている場合は、一旦 $X _ 1, ..., X _ N$ と別々の確率変数が得られる と考えるところからスタートしたのを思い出してください。その時、面倒ではありますが、一旦グラフィカルモデルをひたすらコピーして

f:id:s0sem0y:20200814195031p:plain

と表記したのでした。先程までの流れでは $\mu, sigma$ が共通であることと、$X$ は独立であること(独立同分布から $X$ が生成されること)を仮定し、グラフィカルモデルを簡略したのでした。今回、時系列データの場合は、$X$ が互いに独立であるとは考えられないような場合が多いはずです(そもそも時刻のインデックスとは何なのだろう。すべてのデータは多かれ少なかれ時間のズレがある中でデータが取得されているはずだ。すなわちすべてのデータにタイムスタンプがついており、それは時刻のインデックスを有しているということである。先程の例で我々は、観測データ $X _ i$ が独立同分布から生成されていると考えたことによって、時刻のインデックスは無意味なものであると切り捨てたのである!)。

今回の時刻のインデックスは非常に重要な意味を持つはずであるので、グラフィカルモデルを安易に簡略化するわけにはいきません。一度基本に立ち返り、すべての確率変数を列挙するところから始めます。インデックスを $i$ 改め $t$ として、時刻のインデックスの最大値(最後の時刻)を $T$ と書くことにしましょう。そうして、確率変数を列挙すると $X _ 1, X _ 2, X _ 3, ..., X _ T$ と書くことができます。こうして確率変数を列挙してみて、その確率変数間の関係を線で表現するのでした。今、$X _ t$ という確率変数は1つ前の時刻 $X _ {t - 1}$ と2つ前の $X _ {t - 2}$ に依存することを "仮定" すると

f:id:s0sem0y:20200814200808p:plain

と表記できます。逆に、数式やコトバの説明なしに、下記のようなグラフィカルモデルを見せられた場合には

f:id:s0sem0y:20200814200907p:plain

あ、これは確率変数 $X _ t$ は1つ前の時刻 $X _ {t-1}$ にのみ依存したモデルなのだな、と分かる必要があります。 多くの応用では、この1つ前の時刻に依存するモデルで表現されることが多いです。このようなモデルは1次マルコフモデルと呼びます。一般にn次マルコフモデルは、n個前までの時刻に依存する形式になります。(なので、グラフィカルモデルの1つ目は2次マルコフモデルである)

本当に1つ前よりも更に過去を考慮しなくて良いのだろうか。実は、状態空間モデルというクラスの時系列モデルを構築することで、形式的に1つ前の時刻しか登場しないような表現に変形できる場合が多々ある。状態空間モデルは十分に表現できるクラスが広いため、実用上は1次マルコフモデルの登場がもっとも多い(はず)。このあたりは統計よりもむしろ、制御や信号処理における高階差分方程式の1階連立差分方程式化を参考にすると良い。無論考えうる時系列表現のすべてを状態空間表現に持っていけるわけではないが。

ちなみに、1次マルコフモデルは数式で表現するならば、同時分布の設計の手順を思い出し、

$$ p(X _ 1, X _ 2, ..., X _ t, ..., X _ T) = p(X _ 1)p(X _ 2 \mid X _ 1)...p(X _ t \mid X _ {t - 1})...p(X _ T \mid X _ {T-1}) = p(X _ 1) \prod _ {t = 2} ^ T p(X _ t \mid X _ {t - 1}) $$

と即座に書ける必要があります(もちろん丁寧に読んできていたら、これくらいは簡単に見えているかと思われます)。

このマルコフモデルで表現できる最も単純な例はランダムウォークでしょう。ランダムウォークは1つ前の時刻の値にガウスノイズを乗せて次の値とします。すなわち、

$$ p(X _ t \mid X _ {t - 1}) = \mathcal N(X _ t \mid X _ {t - 1}, \sigma) $$

と書くことができます(平均的には1つ前と同じ値になるのだが、標準偏差 $\sigma$ でバラつくというものである。同時分布には $p(\sigma)$ を追加しなければなりません)。これに伴ってグラフィカルモデルは

f:id:s0sem0y:20200814203405p:plain

そうすると $p(\sigma)$ のハイパーパラメータはどうしよう!とか思えてくるかもしれませんが、これは時系列モデルの例として見せているだけでなのでこれ以上は追う必要はありません。 なぜならランダムウォークは基本的には $X _ t$ の挙動を予測することはできないからです。平均値は、完全に初期の値 $X _ 1$ になるだけでそれ以上の情報はランダムウォークからは得られないのです。すなわち、$p(X _ 1)$ という初期分布が全てとなります(という意味では学習においては楽しい題材ではない。解析対象がランダムウォークだとなったら、基本は尻尾を巻いて逃げるしかないのである)。

最後に

一旦区切ります。

www.hellocybernetics.tech