多観測の尤度計算〜Haskellで確率モデル(続)

  • 昨日の記事メトロポリスヘイスティングをやった
  • その実行例として、二次元点集合に一次線形モデルを設定して回帰係数の事後分布を作った
  • その中でサンプリングされた回帰直線係数の下での条件付き確率を、個別の観察点について計算し、それをたたんで一つの分布関数からの確率にし
  • それをメトロポリスヘイスティングにまわしている
  • そのことをコードに沿って確認する
  • まずはメトロポリスヘイスティング関数から
xs <- Random.sample $ mh 100 100000 (points td linear)
    • xs <- は、<- の右辺がMonadRandomというモナドなので、その中身を取り出している
    • Random.sample は
:t Random.sample
Random.sample :: (MonadRandom m, Distribution d t) => d t -> m t
    • とあるように、Distributionを取って、具体的な乱数を適用した結果tをMonadRandomに包んで返す関数
    • Random.sampleを掛けて初めて、具体的な乱数が使われる。それまでは、いくつの乱数を発生するか、いくつかを捨てるか、どうやって乱数を発生させるか(メトロポリスヘイスティングという仕組みが乱数発生器〜分布〜と見られている)
    • mh 100 100000 (points td linear) はバーンイン100回、記録100000回でメトロポリスヘイスティングサンプリングするという処理で、サンプリングするのは、(points td linear)が作っている分布
  • points td linear。
    • points は、二次元点座標のリストと、linearというDistとを受け取って、Distを作る
    • どんなDistが作りたいかと言えば、(xi, yi)という観察があったときに、平均 axi + b の正規分布においてyiが起きる確率密度を計算し、それをすべての観測点について掛けあわせたようなDist
points :: [Main.Point] -> Dist Param -> Dist Param
points ps d = foldl (flip point) d ps
    • そのような処理をするpointsという関数は上記のように定義されている
    • ps が2次元座標のリスト、dが分布
    • foldl hoge d ps となっているので、点リストが定義されていないときは、d と書かれているlinear という関数が定めるDistそのもの。それは、2つの実数値(ax + bの係数であるa,bの候補値)。それの確率は、と言えば、候補がそれしか無い(pure)なので、確率は1
    • 点リストに一つの点(ps = [p1])があれば、points ps d が作るDistにfoldlは関係なくなり、(flip point) d p1 = point p1 dが返る
    • これは次のpointの定義を読むとわかる
point :: Main.Point -> Dist Param -> Dist Param
point (x,y) = condition (\(a,b) -> Data.Random.pdf (Random.Normal (a*x + b) 1) y)
    • point (x1,y1)がcondition (hoge) であると書かれている。
    • これは、(x1,x2)を使って、「とある正規分布を定義し、そこからy1が得られる確率を返す関数であるところのhogeを受け取った状態の関数であって、さらなる引数Dist Paramの入力受付中」、というものである
    • ここでとある正規分布を定めるために(a,b)を持ったDist Paramも与えないといけないし、それをさらに次のDist Paramとともに活用する予定、という宙ぶらりな状態である
    • 実際、condition=Conditionalであって、これは、以下のように「さっぱりとした書きぶり(何をするのかわからない書きぶり)」のものである
data Dist a where
  Pure :: a -> Dist a
  Bind :: Dist b -> (b -> Dist a) -> Dist a
  Primitive :: RVar a -> Dist a
  Conditional :: (a -> Prob) -> Dist a -> Dist a
  • 代数的データ型を思い出そう
    • Haskellのデータ型は代数的
    • 言い換えると、上記の4つのコンストラクタ(と呼んでよいのだと思うが)で作れるすべての複雑な→構造はすべてDist a である
    • そしてその複雑なものも、prior関数では再帰的にちゃんと対応できるようになっている
    • このcondition=Conditionalが何をするのかは、次の「計算」を教えられて初めて意味を持つ
prior :: Dist a -> Dist (a, Prob)
prior (Conditional c d) = do
  (x,s) <- prior d
  return (x, s * c x)
    • 最後の引数であるd というDistが持っているある確率事象xの確率値 s に、後ろから2番めの引数であるcという関数(事象から確率を返す関数)で、その確率事象xに対応する確率 c xを掛けている
    • 個々の中に、事象xが次々に渡される仕組みが書かれており、この事象がlinearによってサンプリングされた(a,b)のことである
    • したがって、linearがサンプリングした(a,b)について、観測点ごとに axi + b を平均として分散を1とする正規分布がyiを出す確率密度が掛け合わされていくことになる
    • そのぱたぱた
    • これは、点p1の座標(x1,y1)を受け取り、第2引数であるDist(ただしこのDistはParamを持っていて、それがこの場合、線形回帰式の2係数のペア(a,b))の情報を用いることで、平均がa*x+b、分散が1の(一次元)正規分布(Random.Normal (a*x + b) 1)) に基づいて、y1を得るpdf (probability density function)を返す
    • 点リストに2つの点ps=[p1,p2]があるときは、foldl が働き始める
    • 1個目の点p1が作った正規分布とは:『平均がa*x+b、分散が1の(一次元)正規分布(Random.Normal (a*x + b) 1)) に基づいて、y1を得るpdf (probability density function)』
    • この正規分布に第二の点p2を取り出してdと処理したい。それをfoldを使って行うためには、次々に処理するものを第2引数とするような関数が必要
    • pointという関数は第1引数に点を取っているので、このままだと困るから、flipして点を第2引数として取る関数 (flip point)を使う
    • foldにはfoldlとfoldrとがあるが、リストに入っている点を前から取り出すならfoldl
    • さて。最後に点kまででで作り上げたDistにk+1番目の点を取って関数pointを適用したい
    • 点を加えていくと、だんだんにDist d が更新されていく。その更新は、Conditionalコンストラクタによってどんどん複雑な形になっていく(入れ子再帰構造)。ただし、そういっても、Dist Paramであることには変わりなく、そのParamとは、Doubleのペアであって、それは、(a,b)のまま(のようだ)
    • したがって、foldlを繰り返すことで、Dist Paramであるdがどんどん条件付きで長くなる
    • その条件での確率(尤度)を計算するには、priorを使い、そのpriorは長くなったDist Paramでも計算できるように定義されている
    • priorの定義をよく見ると、確率値を持っていないDistを、持っているDistに変える、という関数であることがわかる。長ったらしいcondition関数の入れ子になった関数が、「全観察点の尤度計算用関数」として定義された後で、mhにその関数を渡すと、priorが適用され、proposalという関数の出力からメトロポリスヘイスティングがスタートする。proposalの出力は(a,b)とその尤度。
    • ただし、あいかわらず、この記述も「どうやって計算するか」「どうやって更新するか」を記述しているだけ。実際に乱数を発生させるのは、Random.sampleを使うところ