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

"The Borneo office" by angela7dreams

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

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

さて、この確率分布に従う乱数を発生させたいと思います。今回は一番簡単なメトロポリス法を用いて推移確率を決定します。

Haskellで書きました。自分の実力足らずでいろいろとやばいです。IOモナドのところが何かぎこちないです。実行速度も遅いです。Haskellもっと勉強しないとっ。。

view plainprint?

1. import System.Random  
2. import Text.Printf  
3.   
4. select :: (Double, Double) -> IO (Double, Double)  
5. select (x, y) = do  
6.                     dx <- getStdRandom $ randomR (-5,5)  
7.                     dy <- getStdRandom $ randomR (-5,5)  
8.                     return (x+dx, y+dy)  
9.   
10. prob :: (Double, Double) -> Double  
11. prob (x, y) = 1/(2*pi) * exp(-(x*x+y*y)/2)  
12.   
13. nextStatus :: (Double, Double) -> IO (Double, Double)  
14. nextStatus (x, y) = do  
15.                  (xt, yt) <- select(x, y)  
16.                  let p = min 1  (prob(xt, yt) / prob(x, y))  
17.                  q <- getStdRandom $ randomR (0, 1)  
18.                  if p >= q then return (xt, yt) else return (x, y)  
19.   
20. markovChain :: (Double, Double) -> [IO (Double, Double)]  
21. markovChain (x, y) = iterate (>>= nextStatus) (return (x, y))  
22.   
23. printPoints :: (Double, Double) -> IO()  
24. printPoints (x, y) = printf "%f %f\n" x y  
25.   
26. main = mapM_ (>>= printPoints) . take 3000 . drop 1000 . markovChain $ (0, 0)  

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に近づいていくでしょう。