Ideone.com

fork download

copy

// 『パラレルテンパリングを用いた近似ベイズ計算』 http://b...content-available-to-author-only...n.jp/archives/2901
// 本コードは概ね http://a...content-available-to-author-only...v.org/abs/1108.3423v3 の例題にそって作成していますが，一部異なる点があります。
// 標準出力に混合正規分布 0.45 * N(0, 1^2) + 0.45 * N(0, 0.1^2) + 0.10 * N(5, 1^2) からのランダムサンプルを出力します。
 
open System
 
type Random with
   member this.NextUniform inf sup =
      let u = this.NextDouble ()
      inf + u * (sup - inf)
 
   member this.NextNormal mean variance =
      // Box-Muller 法による正規乱数の生成。
      let u1 = this.NextDouble ()
      let u2 = this.NextDouble ()
      let r = sqrt <| -2.0 * log u1
      let x = r * cos (2.0 * Math.PI * u2)  // 効率は悪いが，面倒なので sin の方は捨てる。
      mean + x * sqrt variance
 
// お題は [-10, 10] だけど [-10, 10) とするが，特に問題はない。
let samplePrior (random : Random) = random.NextUniform -10.0 10.0
 
// サンプルされたパラメーターをもとにデータをシミュレーションする。
let simulateValue (random : Random) theta =
   let r = random.NextDouble ()
   if r < 0.45
   then
      random.NextNormal theta 1.0
   elif r < 0.90
   then
      random.NextNormal theta 0.01
   else
      random.NextNormal (theta - 5.0) 1.0
 
// 観察データを 0.0 として |observed - simulated| を計算する。
let distance x = abs x
 
type AbcResult =
   | Accept of float (* parameter *) * float (* distance *)
   | Reject
 
let rejection random tolerance theta =
   let x = simulateValue random theta
   let d = distance x
   if d < tolerance
   then
      Accept (theta, d)
   else
      Reject
 
let rec rejectionInfinite random tolerance =
   seq {
      yield rejection random tolerance (samplePrior random)
      yield! rejectionInfinite random tolerance
   }
 
// 頻繁にスワップするので配列が速いと思われる。
let swapInPlace i j (array : 'a []) =
   let tmp = array.[i]
   array.[i] <- array.[j]
   array.[j] <- tmp
 
let sampleWithReplacement (random : Random) size (source : 'a []) =
   let length = Array.length source
   let result = Array.zeroCreate size
   for index = 0 to size - 1 do
      let randomIndex = random.Next (length)
      result.[index] <- source.[randomIndex]
   result
 
let random = Random ()
let chainLength = 600000
let currentLength = ref 0
let burnin = 150000
let sampleStep = 500 // ブログでは 100 にしている。
let n = 15
let split length (min, max) =
   let l = length - 1
   let r = Math.Pow (max / min, 1.0 / float l)
   in [| for i in 0 .. l -> min * pown r i |]
let tolerance = split n (0.025, 2.0)
let temperature = split n (1.0, 4.0)
 
// 最初だけ ABC-RS を行う。
let current =
   let chooseAcceptance = Seq.choose (function | Accept (theta, d) -> Some (theta, d) | Reject -> None)
   Array.map (rejectionInfinite random >> chooseAcceptance >> Seq.head) tolerance
incr currentLength
 
let inline q (theta, temperature) = random.NextNormal theta (0.10 * 0.10 * temperature)
let inline isInPriorRange x = -10.0 <= x && x < 10.0
let chainIndexes = [| 0 .. n - 1 |]
while !currentLength < chainLength do
   // ABC-MCMC ステップ
   let theta = Array.zip (Array.map fst current) temperature |> Array.map q
   let values = Array.map (simulateValue random) theta
   // alpha = min { 1, (π(candidate) * q(current -> candidate)) / (π(current) * q(candidate -> current)) } * 1_ρ(values)
   // 1_ρ は indicator。
   // q は平均 0 の正規分布なので比は常に 1 になる。
   // π は一様分布なので範囲外で 0，範囲内で 1 で current は常に範囲内にあることが保障される。
   // したがって
   // * candidate が一様分布の範囲内の時 alpha = 1_ρ(values)
   // * candidate が一様分布の範囲外の時 alpha = 0
   for index = 0 to n - 1 do
      let t = theta.[index]
      let x = values.[index]
      let e = tolerance.[index]
      let d = distance x
      if isInPriorRange t && d < e
      then
         current.[index] <- (t, d)
 
   // 交換ステップ
   for loop = 1 to n do
      let index = sampleWithReplacement random 2 chainIndexes
      Array.sortInPlace index
      let i = Array.min index
      let j = Array.max index
      if i <> j
      then
         let (_, d) = current.[j]
         let e = tolerance.[i]
         if d < e
         then
            swapInPlace i j current
 
   // 出力
   incr currentLength
   if !currentLength > burnin && !currentLength % sampleStep = 0
   then
      printfn "%f" (fst current.[0])

Ly8g44CO44OR44Op44Os44Or44OG44Oz44OR44Oq44Oz44Kw44KS55So44GE44Gf6L+R5Ly844OZ44Kk44K66KiI566X44CPIGh0dHA6Ly9iLi4uY29udGVudC1hdmFpbGFibGUtdG8tYXV0aG9yLW9ubHkuLi5uLmpwL2FyY2hpdmVzLzI5MDEKLy8g5pys44Kz44O844OJ44Gv5qaC44GtIGh0dHA6Ly9hLi4uY29udGVudC1hdmFpbGFibGUtdG8tYXV0aG9yLW9ubHkuLi52Lm9yZy9hYnMvMTEwOC4zNDIzdjMg44Gu5L6L6aGM44Gr44Gd44Gj44Gm5L2c5oiQ44GX44Gm44GE44G+44GZ44GM77yM5LiA6YOo55Ww44Gq44KL54K544GM44GC44KK44G+44GZ44CCCi8vIOaomea6luWHuuWKm+OBq+a3t+WQiOato+imj+WIhuW4gyAwLjQ1ICogTigwLCAxXjIpICsgMC40NSAqIE4oMCwgMC4xXjIpICsgMC4xMCAqIE4oNSwgMV4yKSDjgYvjgonjga7jg6njg7Pjg4Djg6DjgrXjg7Pjg5fjg6vjgpLlh7rlipvjgZfjgb7jgZnjgIIKCm9wZW4gU3lzdGVtCgp0eXBlIFJhbmRvbSB3aXRoCiAgIG1lbWJlciB0aGlzLk5leHRVbmlmb3JtIGluZiBzdXAgPQogICAgICBsZXQgdSA9IHRoaXMuTmV4dERvdWJsZSAoKQogICAgICBpbmYgKyB1ICogKHN1cCAtIGluZikKCiAgIG1lbWJlciB0aGlzLk5leHROb3JtYWwgbWVhbiB2YXJpYW5jZSA9CiAgICAgIC8vIEJveC1NdWxsZXIg5rOV44Gr44KI44KL5q2j6KaP5Lmx5pWw44Gu55Sf5oiQ44CCCiAgICAgIGxldCB1MSA9IHRoaXMuTmV4dERvdWJsZSAoKQogICAgICBsZXQgdTIgPSB0aGlzLk5leHREb3VibGUgKCkKICAgICAgbGV0IHIgPSBzcXJ0IDx8IC0yLjAgKiBsb2cgdTEKICAgICAgbGV0IHggPSByICogY29zICgyLjAgKiBNYXRoLlBJICogdTIpICAvLyDlirnnjofjga/mgqrjgYTjgYzvvIzpnaLlgJLjgarjga7jgacgc2luIOOBruaWueOBr+aNqOOBpuOCi+OAggogICAgICBtZWFuICsgeCAqIHNxcnQgdmFyaWFuY2UKCi8vIOOBiumhjOOBryBbLTEwLCAxMF0g44Gg44GR44GpIFstMTAsIDEwKSDjgajjgZnjgovjgYzvvIznibnjgavllY/poYzjga/jgarjgYTjgIIKbGV0IHNhbXBsZVByaW9yIChyYW5kb20gOiBSYW5kb20pID0gcmFuZG9tLk5leHRVbmlmb3JtIC0xMC4wIDEwLjAKCi8vIOOCteODs+ODl+ODq+OBleOCjOOBn+ODkeODqeODoeODvOOCv+ODvOOCkuOCguOBqOOBq+ODh+ODvOOCv+OCkuOCt+ODn+ODpeODrOODvOOCt+ODp+ODs+OBmeOCi+OAggpsZXQgc2ltdWxhdGVWYWx1ZSAocmFuZG9tIDogUmFuZG9tKSB0aGV0YSA9CiAgIGxldCByID0gcmFuZG9tLk5leHREb3VibGUgKCkKICAgaWYgciA8IDAuNDUKICAgdGhlbgogICAgICByYW5kb20uTmV4dE5vcm1hbCB0aGV0YSAxLjAKICAgZWxpZiByIDwgMC45MAogICB0aGVuCiAgICAgIHJhbmRvbS5OZXh0Tm9ybWFsIHRoZXRhIDAuMDEKICAgZWxzZQogICAgICByYW5kb20uTmV4dE5vcm1hbCAodGhldGEgLSA1LjApIDEuMAoKLy8g6Kaz5a+f44OH44O844K/44KSIDAuMCDjgajjgZfjgaYgfG9ic2VydmVkIC0gc2ltdWxhdGVkfCDjgpLoqIjnrpfjgZnjgovjgIIKbGV0IGRpc3RhbmNlIHggPSBhYnMgeAoKdHlwZSBBYmNSZXN1bHQgPQogICB8IEFjY2VwdCBvZiBmbG9hdCAoKiBwYXJhbWV0ZXIgKikgKiBmbG9hdCAoKiBkaXN0YW5jZSAqKQogICB8IFJlamVjdAoKbGV0IHJlamVjdGlvbiByYW5kb20gdG9sZXJhbmNlIHRoZXRhID0KICAgbGV0IHggPSBzaW11bGF0ZVZhbHVlIHJhbmRvbSB0aGV0YQogICBsZXQgZCA9IGRpc3RhbmNlIHgKICAgaWYgZCA8IHRvbGVyYW5jZQogICB0aGVuCiAgICAgIEFjY2VwdCAodGhldGEsIGQpCiAgIGVsc2UKICAgICAgUmVqZWN0CgpsZXQgcmVjIHJlamVjdGlvbkluZmluaXRlIHJhbmRvbSB0b2xlcmFuY2UgPQogICBzZXEgewogICAgICB5aWVsZCByZWplY3Rpb24gcmFuZG9tIHRvbGVyYW5jZSAoc2FtcGxlUHJpb3IgcmFuZG9tKQogICAgICB5aWVsZCEgcmVqZWN0aW9uSW5maW5pdGUgcmFuZG9tIHRvbGVyYW5jZQogICB9CgovLyDpoLvnuYHjgavjgrnjg6/jg4Pjg5fjgZnjgovjga7jgafphY3liJfjgYzpgJ/jgYTjgajmgJ3jgo/jgozjgovjgIIKbGV0IHN3YXBJblBsYWNlIGkgaiAoYXJyYXkgOiAnYSBbXSkgPQogICBsZXQgdG1wID0gYXJyYXkuW2ldCiAgIGFycmF5LltpXSA8LSBhcnJheS5bal0KICAgYXJyYXkuW2pdIDwtIHRtcAoKbGV0IHNhbXBsZVdpdGhSZXBsYWNlbWVudCAocmFuZG9tIDogUmFuZG9tKSBzaXplIChzb3VyY2UgOiAnYSBbXSkgPQogICBsZXQgbGVuZ3RoID0gQXJyYXkubGVuZ3RoIHNvdXJjZQogICBsZXQgcmVzdWx0ID0gQXJyYXkuemVyb0NyZWF0ZSBzaXplCiAgIGZvciBpbmRleCA9IDAgdG8gc2l6ZSAtIDEgZG8KICAgICAgbGV0IHJhbmRvbUluZGV4ID0gcmFuZG9tLk5leHQgKGxlbmd0aCkKICAgICAgcmVzdWx0LltpbmRleF0gPC0gc291cmNlLltyYW5kb21JbmRleF0KICAgcmVzdWx0CgpsZXQgcmFuZG9tID0gUmFuZG9tICgpCmxldCBjaGFpbkxlbmd0aCA9IDYwMDAwMApsZXQgY3VycmVudExlbmd0aCA9IHJlZiAwCmxldCBidXJuaW4gPSAxNTAwMDAKbGV0IHNhbXBsZVN0ZXAgPSA1MDAgLy8g44OW44Ot44Kw44Gn44GvIDEwMCDjgavjgZfjgabjgYTjgovjgIIKbGV0IG4gPSAxNQpsZXQgc3BsaXQgbGVuZ3RoIChtaW4sIG1heCkgPQogICBsZXQgbCA9IGxlbmd0aCAtIDEKICAgbGV0IHIgPSBNYXRoLlBvdyAobWF4IC8gbWluLCAxLjAgLyBmbG9hdCBsKQogICBpbiBbfCBmb3IgaSBpbiAwIC4uIGwgLT4gbWluICogcG93biByIGkgfF0KbGV0IHRvbGVyYW5jZSA9IHNwbGl0IG4gKDAuMDI1LCAyLjApCmxldCB0ZW1wZXJhdHVyZSA9IHNwbGl0IG4gKDEuMCwgNC4wKQoKLy8g5pyA5Yid44Gg44GRIEFCQy1SUyDjgpLooYzjgYbjgIIKbGV0IGN1cnJlbnQgPQogICBsZXQgY2hvb3NlQWNjZXB0YW5jZSA9IFNlcS5jaG9vc2UgKGZ1bmN0aW9uIHwgQWNjZXB0ICh0aGV0YSwgZCkgLT4gU29tZSAodGhldGEsIGQpIHwgUmVqZWN0IC0+IE5vbmUpCiAgIEFycmF5Lm1hcCAocmVqZWN0aW9uSW5maW5pdGUgcmFuZG9tID4+IGNob29zZUFjY2VwdGFuY2UgPj4gU2VxLmhlYWQpIHRvbGVyYW5jZQppbmNyIGN1cnJlbnRMZW5ndGgKCmxldCBpbmxpbmUgcSAodGhldGEsIHRlbXBlcmF0dXJlKSA9IHJhbmRvbS5OZXh0Tm9ybWFsIHRoZXRhICgwLjEwICogMC4xMCAqIHRlbXBlcmF0dXJlKQpsZXQgaW5saW5lIGlzSW5QcmlvclJhbmdlIHggPSAtMTAuMCA8PSB4ICYmIHggPCAxMC4wCmxldCBjaGFpbkluZGV4ZXMgPSBbfCAwIC4uIG4gLSAxIHxdCndoaWxlICFjdXJyZW50TGVuZ3RoIDwgY2hhaW5MZW5ndGggZG8KICAgLy8gQUJDLU1DTUMg44K544OG44OD44OXCiAgIGxldCB0aGV0YSA9IEFycmF5LnppcCAoQXJyYXkubWFwIGZzdCBjdXJyZW50KSB0ZW1wZXJhdHVyZSB8PiBBcnJheS5tYXAgcQogICBsZXQgdmFsdWVzID0gQXJyYXkubWFwIChzaW11bGF0ZVZhbHVlIHJhbmRvbSkgdGhldGEKICAgLy8gYWxwaGEgPSBtaW4geyAxLCAoz4AoY2FuZGlkYXRlKSAqIHEoY3VycmVudCAtPiBjYW5kaWRhdGUpKSAvICjPgChjdXJyZW50KSAqIHEoY2FuZGlkYXRlIC0+IGN1cnJlbnQpKSB9ICogMV/PgSh2YWx1ZXMpCiAgIC8vIDFfz4Eg44GvIGluZGljYXRvcuOAggogICAvLyBxIOOBr+W5s+WdhyAwIOOBruato+imj+WIhuW4g+OBquOBruOBp+avlOOBr+W4uOOBqyAxIOOBq+OBquOCi+OAggogICAvLyDPgCDjga/kuIDmp5jliIbluIPjgarjga7jgafnr4Tlm7LlpJbjgacgMO+8jOevhOWbsuWGheOBpyAxIOOBpyBjdXJyZW50IOOBr+W4uOOBq+evhOWbsuWGheOBq+OBguOCi+OBk+OBqOOBjOS/nemanOOBleOCjOOCi+OAggogICAvLyDjgZfjgZ/jgYzjgaPjgaYKICAgLy8gKiBjYW5kaWRhdGUg44GM5LiA5qeY5YiG5biD44Gu56+E5Zuy5YaF44Gu5pmCIGFscGhhID0gMV/PgSh2YWx1ZXMpCiAgIC8vICogY2FuZGlkYXRlIOOBjOS4gOanmOWIhuW4g+OBruevhOWbsuWkluOBruaZgiBhbHBoYSA9IDAKICAgZm9yIGluZGV4ID0gMCB0byBuIC0gMSBkbwogICAgICBsZXQgdCA9IHRoZXRhLltpbmRleF0KICAgICAgbGV0IHggPSB2YWx1ZXMuW2luZGV4XQogICAgICBsZXQgZSA9IHRvbGVyYW5jZS5baW5kZXhdCiAgICAgIGxldCBkID0gZGlzdGFuY2UgeAogICAgICBpZiBpc0luUHJpb3JSYW5nZSB0ICYmIGQgPCBlCiAgICAgIHRoZW4KICAgICAgICAgY3VycmVudC5baW5kZXhdIDwtICh0LCBkKQoKICAgLy8g5Lqk5o+b44K544OG44OD44OXCiAgIGZvciBsb29wID0gMSB0byBuIGRvCiAgICAgIGxldCBpbmRleCA9IHNhbXBsZVdpdGhSZXBsYWNlbWVudCByYW5kb20gMiBjaGFpbkluZGV4ZXMKICAgICAgQXJyYXkuc29ydEluUGxhY2UgaW5kZXgKICAgICAgbGV0IGkgPSBBcnJheS5taW4gaW5kZXgKICAgICAgbGV0IGogPSBBcnJheS5tYXggaW5kZXgKICAgICAgaWYgaSA8PiBqCiAgICAgIHRoZW4KICAgICAgICAgbGV0IChfLCBkKSA9IGN1cnJlbnQuW2pdCiAgICAgICAgIGxldCBlID0gdG9sZXJhbmNlLltpXQogICAgICAgICBpZiBkIDwgZQogICAgICAgICB0aGVuCiAgICAgICAgICAgIHN3YXBJblBsYWNlIGkgaiBjdXJyZW50CgogICAvLyDlh7rlipsKICAgaW5jciBjdXJyZW50TGVuZ3RoCiAgIGlmICFjdXJyZW50TGVuZ3RoID4gYnVybmluICYmICFjdXJyZW50TGVuZ3RoICUgc2FtcGxlU3RlcCA9IDAKICAgdGhlbgogICAgICBwcmludGZuICIlZiIgKGZzdCBjdXJyZW50LlswXSk=

Success #stdin #stdout 10.38s 12536KB

stdin

copy

Standard input is empty

stdout

copy

https://ideone.com/KEckA

language:

F# (mono 4.1)

created:

visibility:

public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!

Discover > Sphere Engine API

Discover > IDE Widget

Choose your language