「MCMC のデバッグ」 http://blog.recyclebin.jp/archives/4216

4216.fsx

(*
"Debugging MCMC" http://blog.recyclebin.jp/archives/4216

Usage
------
nuget install -x FsRandom
fsi 4216.fsx
fsi --define:DEBUG 4216.fsx
*)

#if INTERACTIVE
#r "FsRandom/lib/net40/FsRandom.dll"
#endif

open FsRandom
open FsRandom.Statistics

module Seq =
   open Checked
   open LanguagePrimitives

   /// Takes elements by specified number.
   let takeBy n (source:seq<_>) = seq {
      use e = source.GetEnumerator ()
      while e.MoveNext () do
         yield e.Current
         let skip = ref (n - 1)
         while !skip > 0 && e.MoveNext () do
            decr skip
   }

   /// Addtion version of pown.
   let inline private sumn (x : ^a) (n : int) : ^a =
      let rec loop x n result =
         if n <= 0 then
            result
         elif n % 2 = 0 then
            loop (x + x) (n / 2) result
         else
            loop x (n - 1) (result + x)
      loop x n LanguagePrimitives.GenericZero<(^a)>

   /// Calculates variance.
   let inline var (samples : seq<(^a)>) : ^a =
      let count = ref 0
      use enumerator = samples.GetEnumerator ()
      if enumerator.MoveNext () then
         let mutable sum = GenericZero<(^a)>
         let mutable v = GenericZero<(^a)>
         incr count
         sum <- sum + enumerator.Current
         while enumerator.MoveNext () do
            incr count
            let c = enumerator.Current
            sum <- sum + c
            let d = sumn c !count + (-sum)
            v <- v + DivideByInt<(^a)> (d * d) (!count * (!count - 1))
         if !count = 1 then
            invalidArg "samples" "Should contain at least 2 elements."
         else
            DivideByInt<(^a)> v (!count - 1)
      else
         invalidArg "samples" "Empty list."

[<Literal>]
let pi = System.Math.PI
// log-density of normal distribution.
let logNormal mu sd x =
   let u = x - mu
   let v = 2.0 * sd * sd
   -u * u / v - 0.5 * log (v * pi)
 
// The mean value of the prior.
let priorMean = 0.0
// The standard deviation of the prior.
let priorSd = 2.0
 
// Constructs Metropolis-Hastings algorithms using the prior N(0, 2^2) and the proposal N(0, 1^2).
// Small variance is favorable condition for debugging.
let sampler y =
   // Log density of the prior.
   let logPrior mu = logNormal priorMean priorSd mu
   // Log likelihood.
   let logLikelihood mu = Array.sumBy (logNormal mu 1.0)
   // Proposes the next state.
   let propose current = normal (current, 1.0)
   let chain = Seq.markovChain (fun (current, y) -> random {
      let! candidate = propose current
      let alpha =
         // Since the proposal distribution is symmetric, the ratio of the movement is always 1 (0 in log-scale), thus omittable.
         let p mu = logLikelihood mu y + logPrior mu
         exp <| p candidate - p current
      let! u = uniform (0.0, 1.0)
      let next = if u <= alpha then candidate else current
#if DEBUG
      let! y = Array.randomCreate 5 (normal (next, 1.0))
#endif
      return (next, y)
   })
   // Set the initial value of mu as 0 (no matter which one we take).
   chain (0.0, y)
 
// The initial condition of the problem (see the blog).
let y = [|2.9; 4.5; 5.7; 3.9; 4.0|]
let state = createState xorshift (123456789u, 362436069u, 521288629u, 88675123u)
 
let samples =
   sampler y state
   |> Seq.skip 500   // Burn-in period.
   |> Seq.takeBy 10  // To avoid autocorrelation.
   |> Seq.take 1000  // Sample size.
   |> Seq.map fst
   |> Seq.cache
    
printfn "Prior"
printfn "* Mean = %.3f" <| priorMean
printfn "* S.D. = %.3f" <| priorSd
printfn "Sample"
printfn "* Mean = %.3f" <| Seq.average samples
printfn "* S.D. = %.3f" <| sqrt (Seq.var samples)
printfn "Posterior (theory)"
let varTheory = 1.0 / (float (Array.length y) + priorSd ** (-2.0))
let meanTheory = varTheory * Array.sum y + priorMean
printfn "* Mean = %.3f" <| meanTheory
printfn "* S.D. = %.3f" <| sqrt varTheory