Search on the blog

2011年10月23日日曜日

MCMCで二次元正規分布に従う乱数を生成する

"The Borneo office" by angela7dreams

 マルコフ連鎖モンテカルロシミュレーション(MCMC)を用いて、事象(x, y)の発生確率が二次元正規分布



に従う事象をサンプリングしてみます。この分布はグラフに書くと下のようになります。



さて、この確率分布に従う乱数を発生させたいと思います。今回は一番簡単なメトロポリス法を用いて推移確率を決定します。
Haskellで書きました。自分の実力足らずでいろいろとやばいです。IOモナドのところが何かぎこちないです。実行速度も遅いです。Haskellもっと勉強しないとっ。。

import System.Random
import Text.Printf

select :: (Double, Double) -> IO (Double, Double)
select (x, y) = do
dx <- getStdRandom $ randomR (-5,5)
dy <- getStdRandom $ randomR (-5,5)
return (x+dx, y+dy)

prob :: (Double, Double) -> Double
prob (x, y) = 1/(2*pi) * exp(-(x*x+y*y)/2)

nextStatus :: (Double, Double) -> IO (Double, Double)
nextStatus (x, y) = do
(xt, yt) <- select(x, y)
let p = min 1 (prob(xt, yt) / prob(x, y))
q <- getStdRandom $ randomR (0, 1)
if p >= q then return (xt, yt) else return (x, y)

markovChain :: (Double, Double) -> [IO (Double, Double)]
markovChain (x, y) = iterate (>>= nextStatus) (return (x, y))

printPoints :: (Double, Double) -> IO()
printPoints (x, y) = printf "%f %f\n" x y

main = mapM_ (>>= printPoints) . take 3000 . drop 1000 . markovChain $ (0, 0)



それっぽい形になっています。E(x) = -0.00368、E(y) = -0.0121、V(x) = 1.03、V(y) = 1.02でした。試行回数を増やせば、それぞれ0, 0, 1, 1に近づいていくでしょう。

0 件のコメント:

コメントを投稿