Why NaNs?

Main.hs

{-# OPTIONS_GHC -Wall     #-}

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import Control.Monad.Bayes.LogDomain
import Control.Monad.Bayes.Primitive
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Conditional
import Control.Monad.Bayes.Inference

import Numeric.GSL.ODE
import Numeric.LinearAlgebra hiding ( step, Vector )

import qualified Data.Vector as V
import Data.Vector ( Vector )
import Control.Monad

import Control.Monad.Bayes.Sampler

import Debug.Trace

-- number of particles used in PMMH
n_particles :: Int
n_particles = 5

-- model data and constants
h :: Double
h = 0.1

k1', b', d', k2', c' :: Double
k1' = 2.0e2  -- Hare carrying capacity
b'  = 2.0e-2 -- Hare death rate per lynx
d'  = 4.0e-1 -- Lynx death rate
k2' = 2.0e1  -- Lynx carrying capacity
c'  = 4.0e-3 -- Lynx birth rate per hare

-- a :: Double
-- a = 0.5

-- type of simulation state
data S = S {p :: Double, z :: Double, log_alpha :: Double}

parameters :: (MonadDist m, CustomReal m ~ Double) => m (Double, Double)
parameters = do
  mu <- uniform 0 1
  sigma <- uniform 0 5
  return (mu,sigma)

-- initial state of the simulation
initial_state :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> m S
initial_state (mu, sigma) = do
  log_p_init <- normal (log 100) 0.2
  log_z_init <- normal (log  50) 0.1
  log_alpha_init <- normal (log mu) sigma
  return (S (exp log_p_init) (exp log_z_init) log_alpha_init)

-- transition model
transition :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> S -> m S
transition params state = do
  w <- normal 0 (sqrt h)
  -- TODO: ODE solver updates state here
  let a = exp (log_alpha state)
      m = solPp a (p state) (z state)
      newP = m ! 1 ! 0
      newZ = m ! 1 ! 1
      sigma = snd params
      newLog_alpha = -sigma * sigma * h / 2 - sigma * w
  return $ S newP newZ newLog_alpha
    where
      ppOde a k1 b d k2 c _t [pp, zz] =
        [
          a * pp * (1 - pp / k1) - b * pp * zz
        , -d * zz * (1 + zz / k2) + c * pp * zz
        ]
      ppOde _a _k1 _b _d _k2 _c _t vars =
        error $ "ppOde called with: " ++ show (length vars) ++ " variables"

      solPp a x y = odeSolve (ppOde a k1' b' d' k2' c')
                             [x, y]
                             (fromList [0.0, h])

-- full simulation given parameters, returns the parameters
model :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> m (Double, Double)
model params = foldl step (initial_state params) obs >> return params where
  step old_state p_obs = do
    state <- old_state
    new_state <- transition params state
    observe (Continuous (Normal (log (p state)) 0.1)) (log p_obs)
    return new_state

-- full model with particle filter, returns posterior over model parameters
full_model :: (MonadBayes m, CustomReal m ~ Double) => m (Double, Double)
full_model = parameters >>= (collapse . smc (length obs) n_particles . model)

-- PMMH transition kernel
-- monad-bayes does not currently have truncated normals
pmmh_kernel :: (MonadDist m, CustomReal m ~ Double) => [Double] -> m [Double]
pmmh_kernel [mu, sigma] = do
  mu' <- normal mu 1
  sigma' <- normal sigma 1
  return [mu', sigma']
pmmh_kernel xs = error $ "pmmh_kernel called with: " ++ show (length xs) ++ " variables"

-- full PMMH transition step
pmmh_step :: (MonadDist m, CustomReal m ~ Double) => [Double] -> m [Double]
pmmh_step params = do
  params' <- trace ("Params: " ++ show params) $ pmmh_kernel params
  let kernel_density  = unsafeContJointDensity (pmmh_kernel params) params'
  let kernel_density' = unsafeContJointDensity (pmmh_kernel params') params
  pm_density  <- pseudoDensity full_model (map Just params , [])
  pm_density' <- pseudoDensity full_model (map Just params', [])
  let mh_ratio = pm_density' * kernel_density' / (pm_density * kernel_density)
  accept <- bernoulli (min 1 (fromLogDomain mh_ratio))
  trace ("Accept: " ++ show accept) $ return (if accept then params' else params)

iterateNM :: Monad m => Int -> (a -> m a) -> a -> m (Vector a)
iterateNM n f x
  | n == 0 = return $ V.singleton x
  | otherwise = do
      y <- f x
      liftM (y `V.cons`) (iterateNM (n - 1) f y)

main :: IO ()
main = do
  ps1 <- sampleIOfixed $ do ps <- parameters
                            foo <- iterateNM 2 pmmh_step [fst ps, snd ps]
                            return foo
  let mus = V.map (!!0) $ V.drop 1 ps1
  putStrLn $ show (sum mus / 1)
  let sigmas = V.map (!!1) $ V.drop 1 ps1
  putStrLn $ show (sum sigmas / 1)

obs :: [Double]
obs = take 3 $ [
  76.2403679754159,
  65.0098784873532,
  72.7307834964011,
  66.9032651694089,
  69.4511465390588,
  64.2803362618095,
  56.2848299444718,
  55.38783776159,
  58.4038970077208,
  65.1111562800144,
  61.8487358181798,
  54.6268190304816,
  61.7832255452353,
  57.2580992086863,
  51.6506692415061,
  53.513024643873,
  53.2516281763045,
  52.8808601383571,
  50.8172783401389,
  56.4015933848598,
  51.9694880143296,
  53.6405052882993,
  54.575606130029,
  59.3492395433371,
  60.0888983492202,
  44.4110818026195,
  54.4982828776438,
  73.4084712052465,
  54.7644338856397,
  65.9215790688941,
  64.4989883914755,
  80.9772033868803,
  74.4092779374009,
  69.1141635364459,
  67.4281988428452,
  73.2048837536519,
  68.120224957282,
  62.0069646761111,
  80.1892530691497,
  61.1742005270423,
  79.8803349040011,
  83.2200370887438,
  77.5918040799712,
  72.5731170839407,
  77.4547803718811,
  76.6217772509249,
  73.676388095968,
  85.6784703161388,
  104.170224855567,
  89.2220461329041,
  97.8557959956232,
  108.537541058624,
  80.7233691710187,
  101.184198386732,
  106.119103176147,
  100.475879494068,
  125.13032727884,
  91.0701758431548,
  88.8812050725806,
  100.679238971492,
  110.589530407881,
  89.6049957679831,
  103.926588273376,
  100.068112237358,
  115.512346913051,
  113.109941769312,
  114.117083726776,
  126.223426393656,
  116.655177064036,
  119.47938528113,
  119.673062865914,
  122.17342562363,
  115.260394834536,
  118.68195366763,
  124.943355707336,
  151.936975691988,
  114.130007656096,
  148.946645952698,
  144.565611912741,
  118.759397178281,
  126.673416199079,
  142.600114603459,
  124.392431912102,
  123.845615416597,
  131.162727200371,
  127.530085947053,
  149.07114126,
  123.535118463953,
  166.415187417087,
  131.719924167561,
  139.991005535832,
  130.070166508874,
  141.353801701055,
  136.634266000686,
  120.874938386323,
  135.843010192177,
  158.44951183627,
  150.728722863145,
  176.881378505555,
  133.774741754552,
  137.445456445669,
  124.34236836753,
  152.292445573619,
  145.969440861177,
  131.636047513673,
  146.963132155529,
  151.695165885035,
  165.977302905119,
  142.031450539543,
  164.846194851848,
  168.013420195697,
  152.767971695764,
  143.851647250068,
  160.436734291349,
  152.302294176293,
  155.067612585515,
  175.247982500084,
  136.897735270247,
  151.97211044291,
  160.887803973557,
  136.486833382118,
  158.883061549789,
  152.458944513575,
  149.023829176969,
  155.340944696967,
  157.829634728021,
  158.633421934879,
  132.036887061947,
  135.101844477751,
  153.29922033573,
  158.052850728581,
  143.369824868863,
  157.714843025493,
  144.325777788735,
  131.0549741648,
  144.337275322618,
  145.026723249829,
  160.030142868999,
  159.542825278056,
  147.87991947625,
  126.381489671728,
  145.226447360054,
  136.894318063363,
  159.268060107793,
  119.239998434002,
  132.788094593534,
  150.663078585921,
  144.430826746593,
  165.173155950077,
  147.929951578988,
  157.760394555343,
  131.053973697495,
  165.580489628044,
  146.718090811465,
  142.773179403759,
  146.392167581097,
  146.025983883888,
  141.47610191321,
  134.214216068753,
  119.248789990169,
  135.031489547166,
  120.561743419339,
  167.557212051444,
  148.168532426848,
  136.547249833344,
  152.416629324237,
  173.5981255812,
  132.895016695342,
  138.669418114316,
  167.112761599289,
  119.007943156896,
  144.496358594722,
  138.515189246643,
  158.554459930596,
  136.975247626486,
  129.513055665373,
  169.020473191163,
  156.116991826441,
  129.524981775373,
  148.85363969582,
  157.801922552141,
  157.577368682676,
  140.430975061626,
  136.205715759791,
  158.040217033728,
  125.493638253065,
  130.114292183725,
  142.750827896278,
  174.071675002185,
  155.168991141177,
  165.5803346601,
  133.967097475562,
  133.930047497004,
  127.601345216493,
  135.602745924639,
  163.280148579137,
  150.634881342163,
  144.898580660782,
  123.221738509521,
  170.394076047434,
  152.596479493751,
  156.695148025635,
  125.317613954204,
  152.957730717693,
  168.012683439549,
  120.716643641628,
  150.531587838906,
  145.104736014324,
  124.182149449396,
  133.269771329542,
  114.976356844067,
  139.860395302311,
  137.801555771354,
  139.678968034677,
  134.434880549483,
  138.011761267802,
  125.439482936179,
  151.74813237231,
  131.847695399538,
  123.319231908423,
  152.547130543284,
  125.491719298048,
  137.096981763213,
  138.768540718737,
  181.5283483391,
  126.098802734704,
  163.482620601302,
  128.342617368436,
  141.586517467223,
  119.474015267941,
  120.9364132138,
  112.984981940618,
  141.341437153727,
  135.160684269048,
  134.509200686011,
  163.319201985027,
  135.415854723738,
  159.063490681458,
  153.148663733748,
  149.026302196215,
  164.826019093328,
  123.071508271513,
  176.387627783159,
  149.483954577518,
  123.02095153944,
  156.650270662224,
  181.195533023262,
  141.238386891745,
  134.701979169655,
  131.340667101923,
  144.548953780858,
  131.588662675565,
  144.541170469088,
  125.578498701588,
  164.74439785325,
  127.766537075735,
  133.273837701266,
  145.93004226061,
  123.374900614912,
  151.596803755421,
  144.894219796883,
  126.998541138479,
  147.035896225008,
  136.476333130728,
  135.608414577445,
  131.807584253206,
  127.366264864742,
  134.964091100405,
  143.421168142046,
  142.054495699452,
  151.438892001945,
  143.05917211411,
  136.712518442789,
  120.003089311654,
  140.869561010692,
  135.058679734824,
  134.081730533486,
  139.43057733862,
  162.26933796043,
  136.471725650913,
  125.873826665898,
  151.097217418369,
  121.162250703819,
  150.767682408018,
  130.606266701801,
  105.20437547654,
  132.359328972528,
  116.989153980971,
  148.954925443656,
  142.105725642533,
  143.553119435633,
  127.751657026012,
  147.148424004414,
  114.85477950561,
  139.694548684757,
  146.171259764561,
  127.789841217549,
  163.168916396005,
  147.231424950466,
  125.058830151811,
  119.074143576494,
  133.536702795857,
  137.262453768502,
  179.672247639064,
  118.207552485961,
  136.618926147412,
  142.352972731303,
  136.479602357591,
  118.750751452668,
  147.986276828413,
  147.411481380629,
  144.052405584174,
  146.153708105283,
  131.557481669184,
  134.453090884824,
  153.443307349816,
  146.984446479185,
  138.421394583867,
  132.707447826652,
  136.516554357634,
  122.191521371155,
  134.80783216139,
  155.447901778778,
  136.926153057977,
  125.184160305696,
  122.367097582219,
  134.726963096563,
  138.72230909606,
  133.531816345683,
  138.817134160604,
  151.680755115041,
  143.081423812855,
  144.507716078662,
  150.54905673819,
  164.395431529991,
  136.035307291772,
  135.963680724784,
  159.088071380666,
  115.544651209282,
  126.449866575417,
  145.861556307966,
  121.299809587301,
  135.041076468273,
  115.420027815205,
  120.649592931961,
  149.222665704902,
  136.787170312266,
  117.577955349667,
  123.09131466115,
  149.613666057669,
  127.542230989168,
  148.826516120672,
  172.024769211347,
  119.163511852413,
  168.846806388603,
  122.577905606408,
  163.691820942384,
  121.152394606119,
  125.364015734805,
  134.053429227603,
  116.102707413675,
  115.932034184824,
  157.856002757677,
  120.417539311501,
  159.448081983552,
  146.523081800403,
  142.244157886188,
  147.49875516021,
  149.794865803213,
  142.424753997512,
  155.377391656754,
  133.457569330548,
  157.979204500375,
  132.867976837153,
  150.766104208022,
  151.574164752913,
  142.79494077867,
  135.300935952747,
  155.087644429195,
  134.339799297657,
  117.780905975319,
  141.325010750032,
  137.141299531074,
  150.973508888167,
  135.037762357837,
  145.031582826828,
  138.551290507779,
  144.969403511833,
  162.819743874685,
  132.441267184392,
  131.380740568251,
  141.527642238539,
  124.238855213983,
  155.588321253241,
  145.068004565977,
  120.849885116319,
  138.600747823775,
  138.767637532808,
  149.658261434865,
  153.218173715505,
  123.173078726806,
  144.980628125895,
  133.601607881188,
  143.600386923789,
  134.99375013592,
  178.108618865062,
  162.589164120014,
  149.801157966121,
  124.574389244282,
  159.749711094236,
  154.462697945675,
  124.745887246537,
  134.619145288027,
  132.829029331846,
  133.515016847945,
  130.370138175619,
  134.058600191345,
  138.967206000683,
  144.844845114948,
  140.150588785047,
  168.155026956341,
  105.663614240897,
  188.849879696026,
  134.066594679321,
  118.875307519085,
  152.36326566797,
  126.208445705652,
  168.297022523569,
  152.27258877451,
  140.382228072818,
  119.424396673087,
  161.465352774379,
  151.655127807313,
  135.515116635163,
  151.532357393466,
  151.810511351233,
  136.361410128331,
  153.741012514592,
  147.431121676277,
  157.91915964734,
  120.327187813985,
  120.171500909668,
  142.428896431138,
  170.204819149259,
  158.214506571837,
  138.733887337899,
  123.808097192822,
  152.903612727097,
  140.769759838882,
  119.667303148295,
  129.45932949513,
  149.16660912873,
  150.254944400129,
  128.273570440445,
  153.278564595965,
  121.864809541511,
  161.542701435884,
  156.946807512148,
  158.881028441797,
  125.569197718704,
  137.940620707105,
  129.90603859364,
  152.757290961669,
  127.74685675702,
  137.183947553826,
  150.462794723172,
  140.929843406962,
  164.975435125998,
  182.344447097492,
  156.919848441299,
  150.750704137777,
  141.693725748231,
  110.676915031473,
  173.230786634625,
  148.248841760484,
  119.716751121111,
  143.947798672895,
  130.440507876832,
  131.964077776465,
  134.646898631033,
  100.778877180885,
  124.348744604234,
  127.180330688462,
  158.146060814462,
  171.358807156152,
  131.864294544909,
  143.085939246005,
  149.586050733912,
  132.2767374166,
  129.632802049237,
  137.02205312678,
  146.243848392623,
  125.713745961482,
  132.545872818923,
  129.355428299295,
  150.515004650034,
  149.481788573074,
  154.478850748833,
  137.329325663629,
  127.9504469047,
  133.350229974433]