今日 Twitter で見かけたツイート。

せっかくなのでやってみよう。ヒストグラムを出力するようなシミュレーションだし，気軽に乱数計算ができて絵が描ける言語が良いですね。となると，当然 F# でしょう。

NuGet から FSharp.Charting と FsRandom を持ってきます。

nuget install -ExcludeVersion FSharp.Charting
nuget install -ExcludeVersion FsRandom

REPL 環境でやりたいので Fsi を使います。まず必要なライブラリーを参照します。

#r "System.Windows.Forms.DataVisualization.dll"
#r @"FSharp.Charting\lib\net40\FSharp.Charting.dll"
#r @"FsRandom\lib\net40\FsRandom.dll"

open FSharp.Charting
open FsRandom
open FsRandom.Statistics

シミュレーションは単純です。ランダムに異なる 2 要素を選び，一方から他方 1 コインに与えるということを繰り返すだけです。

まずはランダムに異なる 2 要素を選びます。

let takeTwoIndices array = random {
   let size = Array.length array
   let! first = uniformDiscrete (0, size - 1)
   let! second = uniformDiscrete (0, size - 2)
   return first, if first <= second then second + 1 else second
}

そしてコインを与えるという操作です。ある状態に対してコインのやり取りがあって次の状態が決まるのでマルコフ連鎖です。 FsRandom ではマルコフ連鎖も容易に作れます。

let next current = random {
   let! i, j = takeTwoIndices current
   if current.[i] > 0 then
      current.[i] <- current.[i] - 1
      current.[j] <- current.[j] + 1
   return current
}
let simulate = Seq.markovChain next

これだけでシミュレーション部分は完成です。実際にシミュレーションを実行してみましょう。

let seed = 123456789u, 362436069u, 521288629u, 88675123u
let result =
   simulate xorshift seed
   <| Array.create 100 20  // 初期条件
   |> Seq.nth 1000000  // シミュレーションの繰り返し回数

ここではシミュレーション条件を，人数 100，初期コイン 20 枚，シミュレーションの繰り返しかえ数を 100 万回としています。

結果を適当な幅で集計します。

let histogram binWidth distribution =
   let count =
      distribution
      |> Seq.groupBy (fun i -> i / binWidth)
      |> Seq.map (fun (key, s) -> key, Seq.length s)
      |> Map.ofSeq
   let minKey, maxKey =
      Map.toSeq count
      |> Seq.reduce (fun (minKey, maxKey) (key, _) -> min minKey key, max maxKey key)
   [
      for i in minKey .. maxKey ->
         let label = sprintf "[%d,%d)" (i * binWidth) ((i + 1) * binWidth)
         let frequency = Map.tryFind i count |> function | Some (n) -> n | None -> 0
         label, frequency
   ]

手を抜いて階級幅指定して集計するのとラベル付き集計結果を作成するのを 1 つの関数にまとめています。 Map モジュールには reduce 関数がないので一度 seq に戻してから fold を使っています。

最後に集計結果を出力します。集計幅は 4 としています。

let chart = Chart.Column (histogram 4 result)
chart.ShowChart ()

次のような結果を得ます。本当は棒と棒の間の隙間を埋めてちゃんとしたヒストグラムを描きたかったのですが，やり方がよくわかりませんでした。いげ太さんがコメントで教えてくれましたので参考にして図を差し替えました。

最近は R や Python のようにお手軽に計算してグラフ化するということができる言語に需要がありますが， F# だって負けられません。最近は海外の F# 界隈でもそういう動きが活発になっているのを強く感じます。日本もこの波に乗れると良いですね。