乱数のリストを生成する

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

そもそもどうやって乱数をつくるのか
----------------------------
[`mwc-random`](http://hackage.haskell.org/package/mwc-random)というライブラリを使うのが簡単です。実際に使い方を見てみましょう。

```hs: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)
```

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

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

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

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

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

```hs
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
```

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

```hs
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
```

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

```hs: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
```

```bash
$ runhaskell Haskell.hs
3.1196
```

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

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

```hs
$ 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]$ に含まれる点をサンプルできています。次に **サンプリングを繰り返してリストにする関数** を考えましょう

```hs
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`が出来るような型を表しています。これらを組み合せると

```hs
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]`を作ることができます。

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

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

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

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

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

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

```hs: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
```

```bash
$ runhaskell Main.hs
3.112
```

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