乱数のリストを生成する #Haskell

Haskellで乱数のリストを作る方法を書いていきたいと思います。

そもそもどうやって乱数をつくるのか

mwc-randomというライブラリを使うのが簡単です。実際に使い方を見てみましょう。

Main.hs

{-# Language TypeApplications #-}
import System.Random.MWC

main :: IO ()
main = do
  gen <- createSystemRandom
  rand1 <- uniform @Int gen
  rand2 <- uniform @Double gen
  rand3 <- uniformR @Int (1, 10) gen
  print (rand1, rand2, rand3)

$ runhaskell Main.hs
(2380251754465621693,0.7151906692113023,7)

まず createSystemRandom :: IO GenIO を使って乱数のシードとなる値genを作ります。

uniform :: GenIO -> IO a
uniformR :: (a, a) -> GenIO -> IO a

型を制限して書くとuniformとuniformRはこのような型になっていて、それぞれ一様な乱数と指定された範囲の中で一様な乱数を生成することができます。生成する乱数の種類は型によって決まります。GHC の TypeApplications 拡張を利用して指定すると読みやすく書けますね。

乱数のリストを作ってみる

モンテカルロ法で円周率を求めるプログラムを考えてみます。既にサンプリングされた $[0,1]\times[0,1]$ の点列があれば

import Data.List

monteCarloPi :: [(Double, Double)] -> Double
monteCarloPi points =
  let allCount = genericLength points
      isIn (x, y) = x^2 + y^2 < 1
      innerCount = genericLength $ filter isIn points
   in 4 * innerCount / allCount

この様に円周率の近似値が求まるでしょう。本題はどのように点列を生成するかです。まずは愚直に再帰関数を使って作ってみましょう。

samplePoints :: Int -> GenIO -> IO [(Double, Double)]
samplePoints 0 _ = pure []
samplePoints n gen =
  let rand = uniformR @Double (0, 1) gen
      point = (,) <$> rand <*> rand
   in (:) <$> point <*> samplePoints (n-1) gen

以上の関数を組み合わせると

Main.hs

{-# Language TypeApplications #-}
import Data.List
import System.Random.MWC

monteCarloPi :: [(Double, Double)] -> Double
monteCarloPi points =
  let allCount = genericLength points
      isIn (x, y) = x^2 + y^2 < 1
      innerCount = genericLength $ filter isIn points
   in 4 * innerCount / allCount

samplePoints :: Int -> GenIO -> IO [(Double, Double)]
samplePoints 0 _ = pure []
samplePoints n gen =
  let rand = uniformR @Double (0, 1) gen
      point = (,) <$> rand <*> rand
   in (:) <$> point <*> samplePoints (n-1) gen

main :: IO ()
main = do
  gen <- createSystemRandom
  points <- samplePoints 10000 gen
  print $ monteCarloPi points

$ runhaskell Haskell.hs
3.1196

この様になります。
samplePointsを使った方法でも問題はないのですが強いて挙げると再利用性に乏しいところが気になります。理想を言えば 点をサンプルする関数 と サンプリングを繰り返してリストにする関数 に分解して実装できた方がベターでしょう。
そこでまずは 点をサンプルする関数 を実装してみます。

samplePoint :: GenIO -> IO (Double, Double)
samplePoint gen =
  let rand = uniformR @Double (0, 1) gen
   in (,) <$> rand <*> rand

$ ghci
> :l Main.hs
> gen <- createSystemRandom
> samplePoint gen
(0.8032514414885267,0.3392009298592972)
> samplePoint gen
(0.9526691566976283,0.20546775017764862)
> samplePoint gen
(0.3461528929568981,0.4780685125554437)

こんな感じで $[0,1]\times[0,1]$ に含まれる点をサンプルできています。次に サンプリングを繰り返してリストにする関数 を考えましょう

sampling :: Variate a => Int -> (GenIO -> IO a) -> GenIO -> IO [a]
sampling 0 _ _ = pure []
sampling n sample gen = (:) <$> sample gen <*> sampling (n-1) sample gen

Variateはmwc-randomで定義されている型クラスでuniformやuniformRが出来るような型を表しています。これらを組み合せると

main :: IO ()
main = do
  gen <- createSystemRandom
  points <- sampling 10000 samplePoint gen
  print $ monteCarloPi points

このように書けます。これでサンプリング手法を変えることも簡単にできるようになりましたね！

ところでsamplingで結局やりたいことは IO a をたくさん集めて IO [a] にすることです。例えばreplicate :: Int -> a -> [a]を使えばIO aから[IO a]を作ることができます。

$ ghci
> :l Main.hs
> gen <- createSystemRandom
> :t replicate 100 (samplePoint gen)
replicate 100 (samplePoint gen) :: [IO (Double, Double)]

しかし今欲しいのは[IO a]ではなくIO [a]でした。そこで

sequence :: (Traversable t, Monad m) => t (m a) -> m (t a)

を使うとIOと[]の順番をひっくり返すことができます。

> :t sequence $ replicate 100 (samplePoint gen)
sequence $ replicate 100 (samplePoint gen) :: IO [(Double, Double)]

これが欲しかったものでしたね！この様にsequenceとreplicateを組み合わせて自前で実装したsamplingを実現できることがわかりました。最終的にモンテカルロ法で円周率を求めるプログラムは以下のようになりました。

Main.hs

{-# Language TypeApplications #-}
import Data.List
import System.Random.MWC

monteCarloPi :: [(Double, Double)] -> Double
monteCarloPi points =
  let allCount = genericLength points
      isIn (x, y) = x^2 + y^2 < 1
      innerCount = genericLength $ filter isIn points
   in 4 * innerCount / allCount

samplePoint :: GenIO -> IO (Double, Double)
samplePoint gen =
  let rand = uniformR @Double (0, 1) gen
   in (,) <$> rand <*> rand

main :: IO ()
main = do
  gen <- createSystemRandom
  points <- sequence $ replicate 10000 (samplePoint gen)
  print $ monteCarloPi points

$ runhaskell Main.hs
3.112

Haskellではこういった基本的な関数とその組み合わせ方を知れば知るほど、複雑な関数を自分で実装する手間がどんどん省かれてゆくなぁと感じる今日このごろでした。