In [1]:
{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax #-}
import Language.Stochaskell
stochaskell
Stochaskell, version 0.1.0
Copyright (C) 2015-2019 David A Roberts
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistribute it
under certain conditions; see the LICENSE for details.

Using installation directory at 
  /home/jovyan/stochaskell
In [2]:
coalLikelihood :: R->(Z,RVec,RVec)->P RVec
coalLikelihood t (n,s,g) = do
  let rate y = g!(findSortedInsertIndex y s)
      widths = vector
        [ (s!(i+1)) - (s!i)
        | i <- 1...(n-1) ] :: RVec
      w = blockVector [cast (s!1), widths,
        cast (t - (s!n))] :: RVec
      areas = vector [ (w!j) * (g!j)
                     | j <- 1...(n+1) ]
      area = sum' areas
  y <- poissonProcess t rate area
  return y

coalPrior :: R -> P (Z,RVec,RVec)
coalPrior t = do
  let lam = 3; nMax = 30; a = 1; b = 200
  n <- truncated 1 nMax (poisson lam)
  s' <- orderedSample
         (2*n + 1) (uniform 0 t) :: P RVec
  let s = vector [ s'!(2*i) | i <- 1...n ]
  g <- joint vector [ gamma a b
                    | _ <- 1...(n+1) ]
  return (n,s,g)

type Model = (Z,RVec,RVec,RVec)
coal :: R -> P Model
coal t = do
  (n,s,g) <- coalPrior t
  y <- coalLikelihood t (n,s,g)
  return (n,s,g,y)
In [3]:
coalMovePoint :: R -> Model -> P Model
coalMovePoint t (n,s_,g,y) = do
  let s i = if i == 0 then 0
       else if i == (n+1) then t
       else (s_!i)
  j <- uniform 1 n
  u <- uniform 0 1
  let sj' = s (j-1) +
           (s (j+1) - s (j-1)) * u
      s' = s_ `replaceAt` (j, sj')
  return (n,s',g,y)

coalMoveRate :: R -> Model -> P Model
coalMoveRate t (n,s,g,y) = do
  j <- uniform 1 (n+1)
  u <- uniform (-1/2) (1/2)
  let gj = exp (log (g!j) + u)
      g' = g `replaceAt` (j, gj)
  return (n,s,g',y)

coalMoveBirth :: R -> Model -> P Model
coalMoveBirth t (n,s,g,y) = do
  x <- uniform 0 t
  let j = findSortedInsertIndex x s
      s' = s `insertAt` (j, x)
  u <- lognormal 0 1
  let g' = g `insertAt` (j+1, u * (g!j))
  return (n+1,s',g',y)

coalMoveDeath :: R -> Model -> P Model
coalMoveDeath t (n,s,g,y) = do
  j <- uniform 1 n
  let s' = s `deleteAt` j
      g' = g `deleteAt` (j+1)
  return (n-1,s',g',y)
In [4]:
coalMoveProbs :: Z -> (R,R,R,R)
coalMoveProbs n = (e,p,b,d)
  where lam = 3
        nMax = 30
        c = 0.9 * (real lam + 1) / (2 * real lam + 1)
        b = if n == nMax then 0 else if n < (lam - 1) then c else c * real lam / (cast n + 1) :: R
        d = if n == 1    then 0 else if lam < n       then c else c * cast n / real lam :: R
        p = if n == 1 then 0     else (1-b-d)/2 :: R
        e = if n == 1 then 1-b-d else (1-b-d)/2 :: R
In [5]:
coalMove :: R -> Model -> P Model
coalMove t m@(n,_,_,_) = do
  let (e,p,b,d) = coalMoveProbs n
  mixture' [(e, coalMoveRate  t m)
           ,(p, coalMovePoint t m)
           ,(b, coalMoveBirth t m)
           ,(d, coalMoveDeath t m) ]

coalStep :: R -> Model -> IO Model
coalStep t m =
  coal t `rjmc` coalMove t `runStep` m
In [6]:
coalData = [ 1851.203, 1851.632, 1851.969, 1851.975, 1852.314, 1852.347, 1852.358, 1852.385, 1852.977 , 1853.196, 1853.229, 1853.319, 1853.500, 1854.135, 1856.396, 1856.506, 1856.539, 1856.618 , 1857.138, 1857.404, 1857.582, 1858.091, 1858.154, 1858.406, 1858.945, 1860.125, 1860.169 , 1860.591, 1860.851, 1860.919, 1860.971, 1861.185, 1861.738, 1861.836, 1862.138, 1862.893 , 1862.937, 1863.178, 1863.794, 1863.939, 1863.986, 1865.459, 1865.971, 1866.064, 1866.340 , 1866.452, 1866.833, 1866.948, 1866.951, 1867.635, 1867.854, 1867.862, 1868.749, 1868.903 , 1868.988, 1869.251, 1869.442, 1869.554, 1869.809, 1869.875, 1870.124, 1870.515, 1870.559 , 1870.633, 1871.027, 1871.151, 1871.167, 1871.736, 1871.816, 1872.123, 1872.240, 1872.769 , 1873.136, 1874.285, 1874.546, 1874.888, 1874.981, 1875.329, 1875.925, 1875.931, 1875.931 , 1876.966, 1877.064, 1877.105, 1877.190, 1877.779, 1877.809, 1878.184, 1878.195, 1878.236 , 1878.433, 1878.696, 1879.036, 1879.172, 1879.501, 1880.057, 1880.539, 1880.689, 1880.944 , 1881.105, 1881.968, 1882.129, 1882.296, 1882.299, 1882.335, 1882.852, 1883.797, 1883.851 , 1884.073, 1884.856, 1885.168, 1885.464, 1885.979, 1886.617, 1886.693, 1886.754, 1886.921 , 1887.134, 1887.405, 1888.298, 1889.051, 1889.198, 1889.793, 1890.102, 1890.190, 1891.252 , 1891.665, 1892.654, 1893.508, 1894.477, 1895.318, 1896.070, 1896.284, 1896.331, 1899.630 , 1901.393, 1902.671, 1905.056, 1905.188, 1905.524, 1906.773, 1908.136, 1908.270, 1908.629 , 1909.127, 1909.825, 1910.357, 1910.970, 1912.520, 1913.784, 1914.409, 1916.615, 1918.031 , 1922.529, 1922.677, 1923.569, 1927.162, 1928.114, 1930.154, 1930.748, 1931.077, 1931.830 , 1931.884, 1932.065, 1932.864, 1932.875, 1933.883, 1934.723, 1935.643, 1935.695, 1936.596 , 1937.500, 1938.354, 1939.821, 1940.218, 1940.424, 1941.420, 1941.522, 1941.574, 1942.001 , 1942.129, 1942.483, 1946.945, 1947.025, 1947.619, 1947.638, 1947.687, 1951.405, 1957.883 , 1960.489, 1962.220]
In [7]:
coalData' = subtract (head days) <$> days
  where days = [365.2425 * yr | yr <- coalData]
In [8]:
let t = 40907
(n0,s0,g0) <- simulate (coalPrior t)
print (n0,s0,g0)
(4,[5103.668914593937,25700.783608506343,29470.225091303586,36244.890506445045],[4.892522331906208e-3,1.4197347319748954e-2,2.00712577711841e-3,1.3410980131244538e-3,2.615629151408432e-3])
In [9]:
-- automatically generate RJMC transition kernel
let kernel = coal t `rjmc` coalMove t
kernel
\input ->
  let v_0_0 = 1 == input.0_0 :: B
      v_0_1 = 30 == input.0_0 :: B
      v_0_2 = input.0_0 < 2 :: B
      v_0_3 = 1.0 + 3.0 :: R
      v_0_4 = 0.9 * v_0_3 :: R
      v_0_5 = 1.0 + 6.0 :: R
      v_0_6 = v_0_4 / v_0_5 :: R
      v_0_7 = v_0_6 * 3.0 :: R
      v_0_8 = 1.0 + input.0_0 :: R
      v_0_9 = v_0_7 / v_0_8 :: R
      v_0_10 = ifThenElse v_0_2 v_0_6 v_0_9 :: R
      v_0_11 = ifThenElse v_0_1 0.0 v_0_10 :: R
      v_0_12 = 1.0 - v_0_11 :: R
      v_0_13 = 3 < input.0_0 :: B
      v_0_14 = v_0_6 * input.0_0 :: R
      v_0_15 = v_0_14 / 3.0 :: R
      v_0_16 = ifThenElse v_0_13 v_0_6 v_0_15 :: R
      v_0_17 = ifThenElse v_0_0 0.0 v_0_16 :: R
      v_0_18 = v_0_12 - v_0_17 :: R
      v_0_19 = v_0_18 / 2.0 :: R
      v_0_20 = ifThenElse v_0_0 v_0_18 v_0_19 :: R
      v_0_21 = ifThenElse v_0_0 0.0 v_0_19 :: R
      v_0_22 = +s 1 1 1 1 :: Z
      v_0_23 = getExternal x_show_0_0 :: Z
      v_0_24 = 1 + input.0_0 :: Z
      v_0_25 = getExternal x_show_0_1 :: Z
      v_0_26 = getExternal x_show_0_2 :: R
      v_0_27 = getExternal x_show_0_3 :: Z
      v_0_28 = getExternal x_show_0_4 :: R
      v_0_29 = getExternal x_show_0_5 :: R
      v_0_30 = getExternal x_show_0_6 :: R
      v_0_31 = getExternal x_show_0_7 :: Z
      v_0_32 = log input.0_2!x_show_0_1 :: R
      v_0_33 = x_show_0_2 + v_0_32 :: R
      v_0_34 = exp v_0_33 :: R
      v_0_35 = replaceIndex input.0_2 x_show_0_1 v_0_34 :: vector R
      v_0_36 = x_show_0_3 == 1 :: B
      v_0_37 = x_show_0_3 - v_0_24 :: Z
      v_0_38 = v_0_37 == 1 :: B
      v_0_39 = x_show_0_3 - 1 :: Z
      v_0_40 = ifThenElse v_0_38 40907.0 input.0_1!v_0_39 :: R
      v_0_41 = ifThenElse v_0_36 0.0 v_0_40 :: R
      v_0_42 = negate x_show_0_3 :: Z
      v_0_43 = v_0_42 == 1 :: B
      v_0_44 = v_0_24 - x_show_0_3 :: Z
      v_0_45 = v_0_44 == 1 :: B
      v_0_46 = x_show_0_3 + 1 :: Z
      v_0_47 = ifThenElse v_0_45 40907.0 input.0_1!v_0_46 :: R
      v_0_48 = ifThenElse v_0_43 0.0 v_0_47 :: R
      v_0_49 = v_0_48 - v_0_41 :: R
      v_0_50 = v_0_49 * x_show_0_4 :: R
      v_0_51 = v_0_41 + v_0_50 :: R
      v_0_52 = replaceIndex input.0_1 x_show_0_3 v_0_51 :: vector R
      v_0_53 = findSortedInsertIndex x_show_0_5 input.0_1 :: Z
      v_0_54 = insertIndex input.0_1 v_0_53 x_show_0_5 :: vector R
      v_0_55 = v_0_53 + 1 :: Z
      v_0_56 = exp x_show_0_6 :: R
      v_0_57 = v_0_56 * input.0_2!v_0_53 :: R
      v_0_58 = insertIndex input.0_2 v_0_55 v_0_57 :: vector R
      v_0_59 = input.0_0 - 1 :: Z
      v_0_60 = deleteIndex input.0_1 x_show_0_7 :: vector R
      v_0_61 = x_show_0_7 + 1 :: Z
      v_0_62 = deleteIndex input.0_2 v_0_61 :: vector R
      v_0_63 = case x_show_0_0 of
        1 ->
          [input.0_0,input.0_1,v_0_35,input.0_3]
        2 ->
          [input.0_0,v_0_52,input.0_2,input.0_3]
        3 ->
          [v_0_24,v_0_54,v_0_58,input.0_3]
        4 ->
          [v_0_59,v_0_60,v_0_62,input.0_3]
      v_0_64 = poisson_lpdf v_0_63.0_0 3.0 :: R
      v_0_65 = 1 > v_0_63.0_0 :: B
      v_0_66 = 2 * v_0_63.0_0 :: Z
      v_0_67 = v_0_66 + 1 :: Z
      v_0_68 = logFactorial v_0_67 :: R
      v_0_69 = uniform_cdf v_0_63.0_1!1 0.0 40907.0 :: R
      v_0_70 = log v_0_69 :: R
      v_0_71 = v_0_63.0_0 - 1 :: Z
      v_0_72 = 
        [ let v_1_0 = i_1_1 + 1 :: Z
              v_1_1 = uniform_cdf v_0_63.0_1!v_1_0 0.0 40907.0 :: R
              v_1_2 = uniform_cdf v_0_63.0_1!i_1_1 0.0 40907.0 :: R
              v_1_3 = v_1_1 - v_1_2 :: R
              v_1_4 = log v_1_3 :: R
              v_1_5 = 2 * v_1_0 :: Z
              v_1_6 = 2 * i_1_1 :: Z
              v_1_7 = v_1_5 - v_1_6 :: Z
              v_1_8 = v_1_7 - 1 :: Z
              v_1_9 = v_1_4 * v_1_8 :: R
              v_1_10 = logFactorial v_1_8 :: R
              v_1_11 = v_1_9 - v_1_10 :: R
           in v_1_11
        | i_1_1 <- 1...v_0_71 ] :: vector R
      v_0_73 = foldl 0.0 v_0_72 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_74 = uniform_cdf v_0_63.0_1!v_0_63.0_0 0.0 40907.0 :: R
      v_0_75 = 1.0 - v_0_74 :: R
      v_0_76 = log v_0_75 :: R
      v_0_77 = v_0_67 - v_0_66 :: Z
      v_0_78 = v_0_76 * v_0_77 :: R
      v_0_79 = logFactorial v_0_77 :: R
      v_0_80 = v_0_78 - v_0_79 :: R
      v_0_81 = 
        [ let v_1_0 = uniform_lpdf v_0_63.0_1!i_1_1 0.0 40907.0 :: R
           in v_1_0
        | i_1_1 <- 1...v_0_63.0_0 ] :: vector R
      v_0_82 = foldl 0.0 v_0_81 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_83 = +s v_0_68 v_0_70 v_0_73 v_0_80 v_0_82 :: R
      v_0_84 = ifThenElse v_0_65 0.0 v_0_83 :: R
      v_0_85 = vectorSize v_0_63.0_2 :: Z
      v_0_86 = 
        [ i_1_1
        | i_1_1 <- 1...v_0_85 ] :: vector Z
      v_0_87 = foldl 0.0 v_0_86 $ \i_1_12 i_1_11 ->
        let v_1_0 = gamma_lpdf v_0_63.0_2!i_1_11 1.0 200.0 :: R
            v_1_1 = i_1_12 + v_1_0 :: R
         in v_1_1
      v_0_88 = 
        [ let v_1_0 = 2 * i_1_1 :: Z
           in ( let v_0_0 = 2 * i_1_1 :: Z
               in v_0_0 -> let v_0_0 = log input.0_2!x_show_0_1 :: R
                  v_0_1 = x_show_0_2 + v_0_0 :: R
                  v_0_2 = exp v_0_1 :: R
                  v_0_3 = replaceIndex input.0_2 x_show_0_1 v_0_2 :: vector R
                  v_0_4 = x_show_0_3 - 1 :: Z
                  v_0_5 = x_show_0_3 == 1 :: B
                  v_0_6 = 1 + input.0_0 :: Z
                  v_0_7 = x_show_0_3 - v_0_6 :: Z
                  v_0_8 = v_0_7 == 1 :: B
                  v_0_9 = ifThenElse v_0_8 40907.0 input.0_1!v_0_4 :: R
                  v_0_10 = ifThenElse v_0_5 0.0 v_0_9 :: R
                  v_0_11 = x_show_0_3 + 1 :: Z
                  v_0_12 = negate x_show_0_3 :: Z
                  v_0_13 = v_0_12 == 1 :: B
                  v_0_14 = v_0_6 - x_show_0_3 :: Z
                  v_0_15 = v_0_14 == 1 :: B
                  v_0_16 = ifThenElse v_0_15 40907.0 input.0_1!v_0_11 :: R
                  v_0_17 = ifThenElse v_0_13 0.0 v_0_16 :: R
                  v_0_18 = v_0_17 - v_0_10 :: R
                  v_0_19 = v_0_18 * x_show_0_4 :: R
                  v_0_20 = v_0_10 + v_0_19 :: R
                  v_0_21 = replaceIndex input.0_1 x_show_0_3 v_0_20 :: vector R
                  v_0_22 = findSortedInsertIndex x_show_0_5 input.0_1 :: Z
                  v_0_23 = insertIndex input.0_1 v_0_22 x_show_0_5 :: vector R
                  v_0_24 = v_0_22 + 1 :: Z
                  v_0_25 = exp x_show_0_6 :: R
                  v_0_26 = v_0_25 * input.0_2!v_0_22 :: R
                  v_0_27 = insertIndex input.0_2 v_0_24 v_0_26 :: vector R
                  v_0_28 = input.0_0 - 1 :: Z
                  v_0_29 = deleteIndex input.0_1 x_show_0_7 :: vector R
                  v_0_30 = x_show_0_7 + 1 :: Z
                  v_0_31 = deleteIndex input.0_2 v_0_30 :: vector R
                  v_0_32 = case x_show_0_0 of
                    1 ->
                      [input.0_0,input.0_1,v_0_3,input.0_3]
                    2 ->
                      [input.0_0,v_0_21,input.0_2,input.0_3]
                    3 ->
                      [v_0_6,v_0_23,v_0_27,input.0_3]
                    4 ->
                      [v_0_28,v_0_29,v_0_31,input.0_3]
               in v_0_32.0_1!i_1_1 | i_1_1 <- 1...let v_0_0 = log input.0_2!x_show_0_1 :: R
                  v_0_1 = x_show_0_2 + v_0_0 :: R
                  v_0_2 = exp v_0_1 :: R
                  v_0_3 = replaceIndex input.0_2 x_show_0_1 v_0_2 :: vector R
                  v_0_4 = x_show_0_3 - 1 :: Z
                  v_0_5 = x_show_0_3 == 1 :: B
                  v_0_6 = 1 + input.0_0 :: Z
                  v_0_7 = x_show_0_3 - v_0_6 :: Z
                  v_0_8 = v_0_7 == 1 :: B
                  v_0_9 = ifThenElse v_0_8 40907.0 input.0_1!v_0_4 :: R
                  v_0_10 = ifThenElse v_0_5 0.0 v_0_9 :: R
                  v_0_11 = x_show_0_3 + 1 :: Z
                  v_0_12 = negate x_show_0_3 :: Z
                  v_0_13 = v_0_12 == 1 :: B
                  v_0_14 = v_0_6 - x_show_0_3 :: Z
                  v_0_15 = v_0_14 == 1 :: B
                  v_0_16 = ifThenElse v_0_15 40907.0 input.0_1!v_0_11 :: R
                  v_0_17 = ifThenElse v_0_13 0.0 v_0_16 :: R
                  v_0_18 = v_0_17 - v_0_10 :: R
                  v_0_19 = v_0_18 * x_show_0_4 :: R
                  v_0_20 = v_0_10 + v_0_19 :: R
                  v_0_21 = replaceIndex input.0_1 x_show_0_3 v_0_20 :: vector R
                  v_0_22 = findSortedInsertIndex x_show_0_5 input.0_1 :: Z
                  v_0_23 = insertIndex input.0_1 v_0_22 x_show_0_5 :: vector R
                  v_0_24 = v_0_22 + 1 :: Z
                  v_0_25 = exp x_show_0_6 :: R
                  v_0_26 = v_0_25 * input.0_2!v_0_22 :: R
                  v_0_27 = insertIndex input.0_2 v_0_24 v_0_26 :: vector R
                  v_0_28 = input.0_0 - 1 :: Z
                  v_0_29 = deleteIndex input.0_1 x_show_0_7 :: vector R
                  v_0_30 = x_show_0_7 + 1 :: Z
                  v_0_31 = deleteIndex input.0_2 v_0_30 :: vector R
                  v_0_32 = case x_show_0_0 of
                    1 ->
                      [input.0_0,input.0_1,v_0_3,input.0_3]
                    2 ->
                      [input.0_0,v_0_21,input.0_2,input.0_3]
                    3 ->
                      [v_0_6,v_0_23,v_0_27,input.0_3]
                    4 ->
                      [v_0_28,v_0_29,v_0_31,input.0_3]
               in v_0_32.0_0 )!v_1_0
        | i_1_1 <- 1...v_0_63.0_0 ] :: vector R
      v_0_89 = \i_1_1 ->
        let v_1_0 = findSortedInsertIndex i_1_1 v_0_88 :: Z
         in v_0_63.0_2!v_1_0
      v_0_90 = 1 + v_0_63.0_0 :: Z
      v_0_91 = 
        [ let v_1_0 = i_1_1 + 1 :: Z
              v_1_1 = v_0_88!v_1_0 - v_0_88!i_1_1 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_71 ] :: vector R
      v_0_92 = 40907.0 - v_0_88!v_0_63.0_0 :: R
      v_0_93 = +s v_0_71 1 1 :: Z
      v_0_94 = 
        [ let v_1_0 = (block array ([1],[3]) [([1],v_0_88!1),([2],v_0_91),([3],v_0_92)] :: vector R)!i_1_1 * v_0_63.0_2!i_1_1 :: R
           in v_1_0
        | i_1_1 <- 1...v_0_90 ] :: vector R
      v_0_95 = foldl 0.0 v_0_94 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_96 = poissonProcess_lpdf v_0_63.0_3 40907.0 v_0_89 v_0_95 :: R
      v_0_97 = +s v_0_64 v_0_84 v_0_87 v_0_96 :: R
      v_0_98 = poisson_lpdf input.0_0 3.0 :: R
      v_0_99 = 1 > input.0_0 :: B
      v_0_100 = 2 * input.0_0 :: Z
      v_0_101 = v_0_100 + 1 :: Z
      v_0_102 = logFactorial v_0_101 :: R
      v_0_103 = uniform_cdf input.0_1!1 0.0 40907.0 :: R
      v_0_104 = log v_0_103 :: R
      v_0_105 = 
        [ let v_1_0 = i_1_1 + 1 :: Z
              v_1_1 = uniform_cdf input.0_1!v_1_0 0.0 40907.0 :: R
              v_1_2 = uniform_cdf input.0_1!i_1_1 0.0 40907.0 :: R
              v_1_3 = v_1_1 - v_1_2 :: R
              v_1_4 = log v_1_3 :: R
              v_1_5 = 2 * v_1_0 :: Z
              v_1_6 = 2 * i_1_1 :: Z
              v_1_7 = v_1_5 - v_1_6 :: Z
              v_1_8 = v_1_7 - 1 :: Z
              v_1_9 = v_1_4 * v_1_8 :: R
              v_1_10 = logFactorial v_1_8 :: R
              v_1_11 = v_1_9 - v_1_10 :: R
           in v_1_11
        | i_1_1 <- 1...v_0_59 ] :: vector R
      v_0_106 = foldl 0.0 v_0_105 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_107 = uniform_cdf input.0_1!input.0_0 0.0 40907.0 :: R
      v_0_108 = 1.0 - v_0_107 :: R
      v_0_109 = log v_0_108 :: R
      v_0_110 = v_0_101 - v_0_100 :: Z
      v_0_111 = v_0_109 * v_0_110 :: R
      v_0_112 = logFactorial v_0_110 :: R
      v_0_113 = v_0_111 - v_0_112 :: R
      v_0_114 = 
        [ let v_1_0 = uniform_lpdf input.0_1!i_1_1 0.0 40907.0 :: R
           in v_1_0
        | i_1_1 <- 1...input.0_0 ] :: vector R
      v_0_115 = foldl 0.0 v_0_114 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_116 = +s v_0_102 v_0_104 v_0_106 v_0_113 v_0_115 :: R
      v_0_117 = ifThenElse v_0_99 0.0 v_0_116 :: R
      v_0_118 = vectorSize input.0_2 :: Z
      v_0_119 = 
        [ i_1_1
        | i_1_1 <- 1...v_0_118 ] :: vector Z
      v_0_120 = foldl 0.0 v_0_119 $ \i_1_12 i_1_11 ->
        let v_1_0 = gamma_lpdf input.0_2!i_1_11 1.0 200.0 :: R
            v_1_1 = i_1_12 + v_1_0 :: R
         in v_1_1
      v_0_121 = 
        [ let v_1_0 = 2 * i_1_1 :: Z
           in ( let v_0_0 = 2 * i_1_1 :: Z
               in v_0_0 -> input.0_1!i_1_1 | i_1_1 <- 1...input.0_0 )!v_1_0
        | i_1_1 <- 1...input.0_0 ] :: vector R
      v_0_122 = \i_1_1 ->
        let v_1_0 = findSortedInsertIndex i_1_1 v_0_121 :: Z
         in input.0_2!v_1_0
      v_0_123 = 
        [ let v_1_0 = i_1_1 + 1 :: Z
              v_1_1 = v_0_121!v_1_0 - v_0_121!i_1_1 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_59 ] :: vector R
      v_0_124 = 40907.0 - v_0_121!input.0_0 :: R
      v_0_125 = +s v_0_59 1 1 :: Z
      v_0_126 = 
        [ let v_1_0 = (block array ([1],[3]) [([1],v_0_121!1),([2],v_0_123),([3],v_0_124)] :: vector R)!i_1_1 * input.0_2!i_1_1 :: R
           in v_1_0
        | i_1_1 <- 1...v_0_24 ] :: vector R
      v_0_127 = foldl 0.0 v_0_126 $ \i_1_12 i_1_11 ->
        let v_1_0 = i_1_11 + i_1_12 :: R
         in v_1_0
      v_0_128 = poissonProcess_lpdf input.0_3 40907.0 v_0_122 v_0_127 :: R
      v_0_129 = +s v_0_98 v_0_117 v_0_120 v_0_128 :: R
      v_0_130 = v_0_97 - v_0_129 :: R
      v_0_131 = v_0_63.0_1 == input.0_1 :: B
      v_0_132 = v_0_63.0_2 == input.0_2 :: B
      v_0_133 = not v_0_132 :: B
      v_0_134 = v_0_63.0_0 - input.0_0 :: Z
      v_0_135 = v_0_134 == 1 :: B
      v_0_136 = not v_0_135 :: B
      v_0_137 = input.0_0 - v_0_63.0_0 :: Z
      v_0_138 = v_0_137 == 1 :: B
      v_0_139 = not v_0_138 :: B
      v_0_140 = &&s v_0_131 v_0_133 v_0_136 v_0_139 :: B
      v_0_141 = not v_0_140 :: B
      v_0_142 = not v_0_131 :: B
      v_0_143 = &&s v_0_132 v_0_136 v_0_139 v_0_142 :: B
      v_0_144 = not v_0_143 :: B
      v_0_145 = &&s v_0_133 v_0_135 v_0_139 v_0_142 :: B
      v_0_146 = not v_0_145 :: B
      v_0_147 = &&s v_0_133 v_0_136 v_0_138 v_0_142 :: B
      v_0_148 = not v_0_147 :: B
      v_0_149 = &&s v_0_141 v_0_144 v_0_146 v_0_148 :: B
      v_0_150 = 1 == v_0_63.0_0 :: B
      v_0_151 = 30 == v_0_63.0_0 :: B
      v_0_152 = v_0_63.0_0 < 2 :: B
      v_0_153 = 1.0 + v_0_63.0_0 :: R
      v_0_154 = v_0_7 / v_0_153 :: R
      v_0_155 = ifThenElse v_0_152 v_0_6 v_0_154 :: R
      v_0_156 = ifThenElse v_0_151 0.0 v_0_155 :: R
      v_0_157 = 1.0 - v_0_156 :: R
      v_0_158 = 3 < v_0_63.0_0 :: B
      v_0_159 = v_0_6 * v_0_63.0_0 :: R
      v_0_160 = v_0_159 / 3.0 :: R
      v_0_161 = ifThenElse v_0_158 v_0_6 v_0_160 :: R
      v_0_162 = ifThenElse v_0_150 0.0 v_0_161 :: R
      v_0_163 = v_0_157 - v_0_162 :: R
      v_0_164 = v_0_163 / 2.0 :: R
      v_0_165 = ifThenElse v_0_150 v_0_163 v_0_164 :: R
      v_0_166 = ifThenElse v_0_150 0.0 v_0_164 :: R
      v_0_167 = categorical_lpdf 1 (block array ([1],[4]) [([1],v_0_165),([2],v_0_166),([3],v_0_156),([4],v_0_162)] :: Array [(1,v_0_22)] R) :: R
      v_0_168 = categorical_lpdf 2 (block array ([1],[4]) [([1],v_0_165),([2],v_0_166),([3],v_0_156),([4],v_0_162)] :: Array [(1,v_0_22)] R) :: R
      v_0_169 = categorical_lpdf 4 (block array ([1],[4]) [([1],v_0_165),([2],v_0_166),([3],v_0_156),([4],v_0_162)] :: Array [(1,v_0_22)] R) :: R
      v_0_170 = categorical_lpdf 3 (block array ([1],[4]) [([1],v_0_165),([2],v_0_166),([3],v_0_156),([4],v_0_162)] :: Array [(1,v_0_22)] R) :: R
      v_0_171 = ifThenElse v_0_149 0.0 { v_0_140 => v_0_167, v_0_143 => v_0_168, v_0_145 => v_0_169, v_0_147 => v_0_170 } :: R
      v_0_172 = foldr -1 v_0_86 $ \i_1_11 i_1_12 ->
        let v_1_0 = v_0_63.0_2!i_1_11 /= input.0_2!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_173 = discreteUniform_lpdf v_0_172 1 v_0_90 :: R
      v_0_174 = id { v_0_140 => v_0_173, v_0_143 => 0.0, v_0_145 => 0.0, v_0_147 => 0.0 } :: R
      v_0_175 = log input.0_2!v_0_172 :: R
      v_0_176 = log v_0_63.0_2!v_0_172 :: R
      v_0_177 = v_0_175 - v_0_176 :: R
      v_0_178 = uniform_lpdf v_0_177 -0.5 0.5 :: R
      v_0_179 = id { v_0_140 => v_0_178, v_0_143 => 0.0, v_0_145 => 0.0, v_0_147 => 0.0 } :: R
      v_0_180 = vectorSize v_0_63.0_1 :: Z
      v_0_181 = 
        [ i_1_1
        | i_1_1 <- 1...v_0_180 ] :: vector Z
      v_0_182 = foldr -1 v_0_181 $ \i_1_11 i_1_12 ->
        let v_1_0 = v_0_63.0_1!i_1_11 /= input.0_1!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_183 = discreteUniform_lpdf v_0_182 1 v_0_63.0_0 :: R
      v_0_184 = id { v_0_140 => 0.0, v_0_143 => v_0_183, v_0_145 => 0.0, v_0_147 => 0.0 } :: R
      v_0_185 = v_0_182 == 1 :: B
      v_0_186 = v_0_182 - v_0_90 :: Z
      v_0_187 = v_0_186 == 1 :: B
      v_0_188 = v_0_182 - 1 :: Z
      v_0_189 = ifThenElse v_0_187 40907.0 v_0_63.0_1!v_0_188 :: R
      v_0_190 = ifThenElse v_0_185 0.0 v_0_189 :: R
      v_0_191 = input.0_1!v_0_182 - v_0_190 :: R
      v_0_192 = negate v_0_182 :: Z
      v_0_193 = v_0_192 == 1 :: B
      v_0_194 = v_0_90 - v_0_182 :: Z
      v_0_195 = v_0_194 == 1 :: B
      v_0_196 = v_0_182 + 1 :: Z
      v_0_197 = ifThenElse v_0_195 40907.0 v_0_63.0_1!v_0_196 :: R
      v_0_198 = ifThenElse v_0_193 0.0 v_0_197 :: R
      v_0_199 = v_0_198 - v_0_190 :: R
      v_0_200 = v_0_191 / v_0_199 :: R
      v_0_201 = uniform_lpdf v_0_200 0.0 1.0 :: R
      v_0_202 = id { v_0_140 => 0.0, v_0_143 => v_0_201, v_0_145 => 0.0, v_0_147 => 0.0 } :: R
      v_0_203 = v_0_180 + 1 :: Z
      v_0_204 = foldr v_0_203 v_0_181 $ \i_1_11 i_1_12 ->
        let v_1_0 = v_0_63.0_1!i_1_11 /= input.0_1!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_205 = uniform_lpdf input.0_1!v_0_204 0.0 40907.0 :: R
      v_0_206 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_145 => 0.0, v_0_147 => v_0_205 } :: R
      v_0_207 = v_0_85 + 1 :: Z
      v_0_208 = foldr v_0_207 v_0_86 $ \i_1_11 i_1_12 ->
        let v_1_0 = v_0_63.0_2!i_1_11 /= input.0_2!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_209 = findSortedInsertIndex input.0_1!v_0_204 v_0_63.0_1 :: Z
      v_0_210 = input.0_2!v_0_208 / v_0_63.0_2!v_0_209 :: R
      v_0_211 = log v_0_210 :: R
      v_0_212 = normal_lpdf v_0_211 0.0 1.0 :: R
      v_0_213 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_145 => 0.0, v_0_147 => v_0_212 } :: R
      v_0_214 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_145 => 0.0, v_0_147 => v_0_213 } :: R
      v_0_215 = vectorSize input.0_1 :: Z
      v_0_216 = v_0_215 + 1 :: Z
      v_0_217 = 
        [ i_1_1
        | i_1_1 <- 1...v_0_215 ] :: vector Z
      v_0_218 = foldr v_0_216 v_0_217 $ \i_1_11 i_1_12 ->
        let v_1_0 = input.0_1!i_1_11 /= v_0_63.0_1!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_219 = discreteUniform_lpdf v_0_218 1 v_0_63.0_0 :: R
      v_0_220 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_145 => v_0_219, v_0_147 => 0.0 } :: R
      v_0_221 = +s v_0_171 v_0_174 v_0_179 v_0_184 v_0_202 v_0_206 v_0_214 v_0_220 :: R
      v_0_222 = categorical_lpdf 1 (block array ([1],[4]) [([1],v_0_20),([2],v_0_21),([3],v_0_11),([4],v_0_17)] :: Array [(1,v_0_22)] R) :: R
      v_0_223 = categorical_lpdf 2 (block array ([1],[4]) [([1],v_0_20),([2],v_0_21),([3],v_0_11),([4],v_0_17)] :: Array [(1,v_0_22)] R) :: R
      v_0_224 = categorical_lpdf 4 (block array ([1],[4]) [([1],v_0_20),([2],v_0_21),([3],v_0_11),([4],v_0_17)] :: Array [(1,v_0_22)] R) :: R
      v_0_225 = categorical_lpdf 3 (block array ([1],[4]) [([1],v_0_20),([2],v_0_21),([3],v_0_11),([4],v_0_17)] :: Array [(1,v_0_22)] R) :: R
      v_0_226 = ifThenElse v_0_149 0.0 { v_0_140 => v_0_222, v_0_143 => v_0_223, v_0_147 => v_0_224, v_0_145 => v_0_225 } :: R
      v_0_227 = foldr -1 v_0_119 $ \i_1_11 i_1_12 ->
        let v_1_0 = input.0_2!i_1_11 /= v_0_63.0_2!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_228 = discreteUniform_lpdf v_0_227 1 v_0_24 :: R
      v_0_229 = id { v_0_140 => v_0_228, v_0_143 => 0.0, v_0_147 => 0.0, v_0_145 => 0.0 } :: R
      v_0_230 = log v_0_63.0_2!v_0_227 :: R
      v_0_231 = log input.0_2!v_0_227 :: R
      v_0_232 = v_0_230 - v_0_231 :: R
      v_0_233 = uniform_lpdf v_0_232 -0.5 0.5 :: R
      v_0_234 = id { v_0_140 => v_0_233, v_0_143 => 0.0, v_0_147 => 0.0, v_0_145 => 0.0 } :: R
      v_0_235 = foldr -1 v_0_217 $ \i_1_11 i_1_12 ->
        let v_1_0 = input.0_1!i_1_11 /= v_0_63.0_1!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_236 = discreteUniform_lpdf v_0_235 1 input.0_0 :: R
      v_0_237 = id { v_0_140 => 0.0, v_0_143 => v_0_236, v_0_147 => 0.0, v_0_145 => 0.0 } :: R
      v_0_238 = v_0_235 == 1 :: B
      v_0_239 = v_0_235 - v_0_24 :: Z
      v_0_240 = v_0_239 == 1 :: B
      v_0_241 = v_0_235 - 1 :: Z
      v_0_242 = ifThenElse v_0_240 40907.0 input.0_1!v_0_241 :: R
      v_0_243 = ifThenElse v_0_238 0.0 v_0_242 :: R
      v_0_244 = v_0_63.0_1!v_0_235 - v_0_243 :: R
      v_0_245 = negate v_0_235 :: Z
      v_0_246 = v_0_245 == 1 :: B
      v_0_247 = v_0_24 - v_0_235 :: Z
      v_0_248 = v_0_247 == 1 :: B
      v_0_249 = v_0_235 + 1 :: Z
      v_0_250 = ifThenElse v_0_248 40907.0 input.0_1!v_0_249 :: R
      v_0_251 = ifThenElse v_0_246 0.0 v_0_250 :: R
      v_0_252 = v_0_251 - v_0_243 :: R
      v_0_253 = v_0_244 / v_0_252 :: R
      v_0_254 = uniform_lpdf v_0_253 0.0 1.0 :: R
      v_0_255 = id { v_0_140 => 0.0, v_0_143 => v_0_254, v_0_147 => 0.0, v_0_145 => 0.0 } :: R
      v_0_256 = uniform_lpdf v_0_63.0_1!v_0_218 0.0 40907.0 :: R
      v_0_257 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_147 => 0.0, v_0_145 => v_0_256 } :: R
      v_0_258 = v_0_118 + 1 :: Z
      v_0_259 = foldr v_0_258 v_0_119 $ \i_1_11 i_1_12 ->
        let v_1_0 = input.0_2!i_1_11 /= v_0_63.0_2!i_1_11 :: B
            v_1_1 = ifThenElse v_1_0 i_1_11 i_1_12 :: Z
         in v_1_1
      v_0_260 = findSortedInsertIndex v_0_63.0_1!v_0_218 input.0_1 :: Z
      v_0_261 = v_0_63.0_2!v_0_259 / input.0_2!v_0_260 :: R
      v_0_262 = log v_0_261 :: R
      v_0_263 = normal_lpdf v_0_262 0.0 1.0 :: R
      v_0_264 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_147 => 0.0, v_0_145 => v_0_263 } :: R
      v_0_265 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_147 => 0.0, v_0_145 => v_0_264 } :: R
      v_0_266 = discreteUniform_lpdf v_0_204 1 input.0_0 :: R
      v_0_267 = id { v_0_140 => 0.0, v_0_143 => 0.0, v_0_147 => v_0_266, v_0_145 => 0.0 } :: R
      v_0_268 = +s v_0_226 v_0_229 v_0_234 v_0_237 v_0_255 v_0_257 v_0_265 v_0_267 :: R
      v_0_269 = v_0_221 - v_0_268 :: R
      v_0_270 = eye v_0_215 :: matrix R
      v_0_271 = ifThenElse { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => v_0_240 } 40907.0 { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => input.0_1!v_0_241 } :: R
      v_0_272 = ifThenElse { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => v_0_238 } 0.0 v_0_271 :: R
      v_0_273 = ifThenElse { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => v_0_248 } 40907.0 { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => input.0_1!v_0_249 } :: R
      v_0_274 = ifThenElse { { v_0_140 => 0, v_0_143 => 1, v_0_147 => 0, v_0_145 => 0 } => v_0_246 } 0.0 v_0_273 :: R
      v_0_275 = v_0_274 - v_0_272 :: R
      v_0_276 = v_0_275 * v_0_253 :: R
      v_0_277 = v_0_272 + v_0_276 :: R
      v_0_278 = replaceIndex input.0_1 v_0_235 v_0_277 :: vector R
      v_0_279 = vectorSize v_0_278 :: Z
      v_0_280 = deleteIndex input.0_1 v_0_204 :: vector R
      v_0_281 = vectorSize v_0_280 :: Z
      v_0_282 = insertIndex input.0_1 v_0_260 v_0_63.0_1!v_0_218 :: vector R
      v_0_283 = vectorSize v_0_282 :: Z
      v_0_284 = zeros v_0_215 v_0_118 :: matrix R
      v_0_285 = vectorSize input.0_3 :: Z
      v_0_286 = zeros v_0_215 v_0_285 :: matrix R
      v_0_287 = v_0_231 + v_0_232 :: R
      v_0_288 = exp v_0_287 :: R
      v_0_289 = replaceIndex input.0_2 v_0_227 v_0_288 :: vector R
      v_0_290 = vectorSize v_0_289 :: Z
      v_0_291 = v_0_204 + 1 :: Z
      v_0_292 = deleteIndex input.0_2 v_0_291 :: vector R
      v_0_293 = vectorSize v_0_292 :: Z
      v_0_294 = v_0_260 + 1 :: Z
      v_0_295 = v_0_261 * input.0_2!v_0_260 :: R
      v_0_296 = insertIndex input.0_2 v_0_294 v_0_295 :: vector R
      v_0_297 = vectorSize v_0_296 :: Z
      v_0_298 = zeros v_0_290 v_0_215 :: matrix R
      v_0_299 = eye v_0_118 :: matrix R
      v_0_300 = zeros 1 v_0_118 :: matrix R
      v_0_301 = 
        [ let v_1_0 = i_1_1 == v_0_227 :: B
              v_1_1 = ifThenElse { v_0_140 => v_1_0 } 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_118 ] :: vector R
      v_0_302 = asRow v_0_301 :: matrix R
      v_0_303 = v_0_302 / input.0_2!v_0_227 :: matrix R
      v_0_304 = v_0_300 + v_0_303 :: matrix R
      v_0_305 = v_0_288 * v_0_304 :: matrix R
      v_0_306 = replaceIndex v_0_299 v_0_227 v_0_305!1 :: matrix R
      v_0_307 = zeros v_0_290 v_0_285 :: matrix R
      v_0_308 = zeros v_0_285 v_0_215 :: matrix R
      v_0_309 = zeros v_0_285 v_0_118 :: matrix R
      v_0_310 = eye v_0_285 :: matrix R
      v_0_311 = +s v_0_118 v_0_215 v_0_285 :: Z
      v_0_312 = zeros v_0_215 1 :: matrix R
      v_0_313 = zeros v_0_118 1 :: matrix R
      v_0_314 = v_0_288 * [[1]] :: matrix R
      v_0_315 = replaceIndex v_0_313 v_0_227 v_0_314!1 :: matrix R
      v_0_316 = zeros v_0_285 1 :: matrix R
      v_0_317 = zeros 1 v_0_215 :: matrix R
      v_0_318 = matrixCols v_0_270 :: Z
      v_0_319 = zeros 1 v_0_318 :: matrix R
      v_0_320 = +s v_0_317 v_0_317 v_0_317 v_0_319 :: matrix R
      v_0_321 = 
        [ let v_1_0 = i_1_1 == v_0_172 :: B
              v_1_1 = ifThenElse v_1_0 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_85 ] :: vector R
      v_0_322 = asRow v_0_321 :: matrix R
      v_0_323 = v_0_322 / v_0_63.0_2!v_0_172 :: matrix R
      v_0_324 = negate v_0_323 :: matrix R
      v_0_325 = v_0_324 <> v_0_306 :: matrix R
      v_0_326 = 
        [ let v_1_0 = i_1_1 == v_0_172 :: B
              v_1_1 = ifThenElse v_1_0 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_118 ] :: vector R
      v_0_327 = asRow v_0_326 :: matrix R
      v_0_328 = v_0_327 / input.0_2!v_0_172 :: matrix R
      v_0_329 = +s v_0_300 v_0_300 v_0_325 v_0_328 :: matrix R
      v_0_330 = zeros 1 v_0_285 :: matrix R
      v_0_331 = matrixCols v_0_310 :: Z
      v_0_332 = zeros 1 v_0_331 :: matrix R
      v_0_333 = +s v_0_330 v_0_330 v_0_330 v_0_332 :: matrix R
      v_0_334 = matrixCols v_0_312 :: Z
      v_0_335 = zeros 1 v_0_334 :: matrix R
      v_0_336 = v_0_324 <> v_0_315 :: matrix R
      v_0_337 = matrixCols v_0_316 :: Z
      v_0_338 = zeros 1 v_0_337 :: matrix R
      v_0_339 = +s v_0_335 v_0_336 v_0_338 :: matrix R
      v_0_340 = +s v_0_118 v_0_215 v_0_285 1 :: Z
      v_0_341 = 
        [ let v_1_0 = v_0_235 - i_1_1 :: Z
              v_1_1 = v_1_0 == 1 :: B
              v_1_2 = ifThenElse { v_0_143 => v_1_1 } 1.0 0.0 :: R
           in v_1_2
        | i_1_1 <- 1...v_0_215 ] :: vector R
      v_0_342 = asRow v_0_341 :: matrix R
      v_0_343 = ifThenElse { v_0_143 => v_0_240 } v_0_317 v_0_342 :: matrix R
      v_0_344 = ifThenElse { v_0_143 => v_0_238 } v_0_317 v_0_343 :: matrix R
      v_0_345 = 
        [ let v_1_0 = i_1_1 - v_0_235 :: Z
              v_1_1 = v_1_0 == 1 :: B
              v_1_2 = ifThenElse { v_0_143 => v_1_1 } 1.0 0.0 :: R
           in v_1_2
        | i_1_1 <- 1...v_0_215 ] :: vector R
      v_0_346 = asRow v_0_345 :: matrix R
      v_0_347 = ifThenElse { v_0_143 => v_0_248 } v_0_317 v_0_346 :: matrix R
      v_0_348 = ifThenElse { v_0_143 => v_0_246 } v_0_317 v_0_347 :: matrix R
      v_0_349 = v_0_348 - v_0_344 :: matrix R
      v_0_350 = v_0_349 * v_0_253 :: matrix R
      v_0_351 = +s v_0_317 v_0_344 v_0_350 :: matrix R
      v_0_352 = replaceIndex v_0_270 v_0_235 v_0_351!1 :: matrix R
      v_0_353 = ifThenElse { v_0_143 => v_0_240 } 40907.0 { v_0_143 => input.0_1!v_0_241 } :: R
      v_0_354 = ifThenElse { v_0_143 => v_0_238 } 0.0 v_0_353 :: R
      v_0_355 = ifThenElse { v_0_143 => v_0_248 } 40907.0 { v_0_143 => input.0_1!v_0_249 } :: R
      v_0_356 = ifThenElse { v_0_143 => v_0_246 } 0.0 v_0_355 :: R
      v_0_357 = v_0_356 - v_0_354 :: R
      v_0_358 = v_0_357 * v_0_253 :: R
      v_0_359 = v_0_354 + v_0_358 :: R
      v_0_360 = replaceIndex input.0_1 v_0_235 v_0_359 :: vector R
      v_0_361 = vectorSize v_0_360 :: Z
      v_0_362 = zeros v_0_361 v_0_118 :: matrix R
      v_0_363 = zeros v_0_361 v_0_285 :: matrix R
      v_0_364 = zeros v_0_118 v_0_215 :: matrix R
      v_0_365 = zeros v_0_118 v_0_285 :: matrix R
      v_0_366 = v_0_357 * [[1]] :: matrix R
      v_0_367 = replaceIndex v_0_312 v_0_235 v_0_366!1 :: matrix R
      v_0_368 = 
        [ let v_1_0 = i_1_1 == v_0_182 :: B
              v_1_1 = ifThenElse v_1_0 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_215 ] :: vector R
      v_0_369 = asRow v_0_368 :: matrix R
      v_0_370 = v_0_369 * v_0_199 :: matrix R
      v_0_371 = v_0_199 ** 2 :: R
      v_0_372 = v_0_370 / v_0_371 :: matrix R
      v_0_373 = zeros 1 v_0_180 :: matrix R
      v_0_374 = 
        [ let v_1_0 = v_0_182 - i_1_1 :: Z
              v_1_1 = v_1_0 == 1 :: B
              v_1_2 = ifThenElse v_1_1 1.0 0.0 :: R
           in v_1_2
        | i_1_1 <- 1...v_0_180 ] :: vector R
      v_0_375 = asRow v_0_374 :: matrix R
      v_0_376 = ifThenElse v_0_187 v_0_373 v_0_375 :: matrix R
      v_0_377 = ifThenElse v_0_185 v_0_373 v_0_376 :: matrix R
      v_0_378 = negate v_0_377 :: matrix R
      v_0_379 = v_0_378 * v_0_199 :: matrix R
      v_0_380 = 
        [ let v_1_0 = i_1_1 - v_0_182 :: Z
              v_1_1 = v_1_0 == 1 :: B
              v_1_2 = ifThenElse v_1_1 1.0 0.0 :: R
           in v_1_2
        | i_1_1 <- 1...v_0_180 ] :: vector R
      v_0_381 = asRow v_0_380 :: matrix R
      v_0_382 = ifThenElse v_0_195 v_0_373 v_0_381 :: matrix R
      v_0_383 = ifThenElse v_0_193 v_0_373 v_0_382 :: matrix R
      v_0_384 = v_0_383 - v_0_377 :: matrix R
      v_0_385 = v_0_191 * v_0_384 :: matrix R
      v_0_386 = v_0_379 - v_0_385 :: matrix R
      v_0_387 = v_0_386 / v_0_371 :: matrix R
      v_0_388 = v_0_387 <> v_0_352 :: matrix R
      v_0_389 = +s v_0_317 v_0_317 v_0_372 v_0_388 :: matrix R
      v_0_390 = matrixCols v_0_299 :: Z
      v_0_391 = zeros 1 v_0_390 :: matrix R
      v_0_392 = +s v_0_300 v_0_300 v_0_300 v_0_391 :: matrix R
      v_0_393 = matrixCols v_0_313 :: Z
      v_0_394 = zeros 1 v_0_393 :: matrix R
      v_0_395 = v_0_387 <> v_0_367 :: matrix R
      v_0_396 = +s v_0_338 v_0_394 v_0_395 :: matrix R
      v_0_397 = v_0_215 - 1 :: Z
      v_0_398 = deleteIndex v_0_270 v_0_204 :: matrix R
      v_0_399 = zeros v_0_281 v_0_118 :: matrix R
      v_0_400 = zeros v_0_281 v_0_285 :: matrix R
      v_0_401 = zeros v_0_293 v_0_215 :: matrix R
      v_0_402 = v_0_118 - 1 :: Z
      v_0_403 = deleteIndex v_0_299 v_0_291 :: matrix R
      v_0_404 = zeros v_0_293 v_0_285 :: matrix R
      v_0_405 = matrixCols v_0_398 :: Z
      v_0_406 = zeros 1 v_0_405 :: matrix R
      v_0_407 = 
        [ let v_1_0 = i_1_1 == v_0_204 :: B
              v_1_1 = ifThenElse v_1_0 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_215 ] :: vector R
      v_0_408 = asRow v_0_407 :: matrix R
      v_0_409 = +s v_0_317 v_0_317 v_0_406 v_0_408 :: matrix R
      v_0_410 = matrixCols v_0_403 :: Z
      v_0_411 = zeros 1 v_0_410 :: matrix R
      v_0_412 = +s v_0_300 v_0_300 v_0_300 v_0_411 :: matrix R
      v_0_413 = +s v_0_317 v_0_317 v_0_317 v_0_406 :: matrix R
      v_0_414 = 
        [ let v_1_0 = i_1_1 == v_0_208 :: B
              v_1_1 = ifThenElse v_1_0 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_118 ] :: vector R
      v_0_415 = asRow v_0_414 :: matrix R
      v_0_416 = v_0_415 * v_0_63.0_2!v_0_209 :: matrix R
      v_0_417 = v_0_63.0_2!v_0_209 ** 2 :: R
      v_0_418 = v_0_416 / v_0_417 :: matrix R
      v_0_419 = v_0_418 / v_0_210 :: matrix R
      v_0_420 = 
        [ let v_1_0 = i_1_1 == v_0_209 :: B
              v_1_1 = ifThenElse { v_0_147 => v_1_0 } 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_85 ] :: vector R
      v_0_421 = asRow v_0_420 :: matrix R
      v_0_422 = input.0_2!v_0_208 * v_0_421 :: matrix R
      v_0_423 = negate v_0_422 :: matrix R
      v_0_424 = v_0_423 / v_0_417 :: matrix R
      v_0_425 = v_0_424 / v_0_210 :: matrix R
      v_0_426 = v_0_425 <> v_0_403 :: matrix R
      v_0_427 = +s v_0_300 v_0_300 v_0_419 v_0_426 :: matrix R
      v_0_428 = insertIndex v_0_270 v_0_260 v_0_317!1 :: matrix R
      v_0_429 = zeros v_0_283 v_0_118 :: matrix R
      v_0_430 = zeros v_0_283 v_0_285 :: matrix R
      v_0_431 = insertIndex v_0_312 v_0_260 [[1]]!1 :: matrix R
      v_0_432 = insertIndex v_0_312 v_0_260 [[0]]!1 :: matrix R
      v_0_433 = zeros v_0_297 v_0_215 :: matrix R
      v_0_434 = 
        [ let v_1_0 = i_1_1 == v_0_260 :: B
              v_1_1 = ifThenElse { v_0_145 => v_1_0 } 1.0 0.0 :: R
           in v_1_1
        | i_1_1 <- 1...v_0_118 ] :: vector R
      v_0_435 = asRow v_0_434 :: matrix R
      v_0_436 = v_0_261 * v_0_435 :: matrix R
      v_0_437 = v_0_300 + v_0_436 :: matrix R
      v_0_438 = insertIndex v_0_299 v_0_294 v_0_437!1 :: matrix R
      v_0_439 = zeros v_0_297 v_0_285 :: matrix R
      v_0_440 = insertIndex v_0_313 v_0_294 [[0]]!1 :: matrix R
      v_0_441 = v_0_261 * [[1]] :: matrix R
      v_0_442 = v_0_441 * input.0_2!v_0_260 :: matrix R
      v_0_443 = insertIndex v_0_313 v_0_294 v_0_442!1 :: matrix R
      v_0_444 = +s v_0_118 v_0_215 v_0_285 1 1 :: Z
      v_0_445 = +s v_0_118 v_0_215 v_0_285 1 1 1 1 :: Z
      v_0_446 = log_det { v_0_140 => block[[block[[v_0_270,v_0_284,v_0_286],[v_0_298,v_0_306,v_0_307],[v_0_308,v_0_309,v_0_310]],block[[v_0_312],[v_0_315],[v_0_316]]],[block[[v_0_320,v_0_329,v_0_333]],v_0_339]], v_0_143 => block[[block[[v_0_352,v_0_362,v_0_363],[v_0_364,v_0_299,v_0_365],[v_0_308,v_0_309,v_0_310]],block[[v_0_367],[v_0_313],[v_0_316]]],[block[[v_0_389,v_0_392,v_0_333]],v_0_396]], v_0_147 => block[[v_0_398,v_0_399,v_0_400],[v_0_401,v_0_403,v_0_404],[v_0_308,v_0_309,v_0_310],[v_0_409,v_0_412,v_0_333],[v_0_413,v_0_427,v_0_333]], v_0_145 => block[[v_0_428,v_0_429,v_0_430,v_0_431,v_0_432],[v_0_433,v_0_438,v_0_439,v_0_440,v_0_443],[v_0_308,v_0_309,v_0_310,v_0_316,v_0_316]] } :: R
      v_0_447 = v_0_269 + v_0_446 :: R
      v_0_448 = +s v_0_130 v_0_269 v_0_446 :: R
      v_0_449 = exp v_0_448 :: R
      v_0_450 = min 1.0 v_0_449 :: R
      v_0_451 = getExternal x_show_0_8 :: B
      v_0_452 = ifThenElse x_show_0_8 v_0_63.0_0 input.0_0 :: Z
      v_0_453 = ifThenElse x_show_0_8 v_0_63.0_1 input.0_1 :: vector R
      v_0_454 = ifThenElse x_show_0_8 v_0_63.0_2 input.0_2 :: vector R
      v_0_455 = ifThenElse x_show_0_8 v_0_63.0_3 input.0_3 :: vector R
   in do x_show_0_0 <- categorical (block array ([1],[4]) [([1],v_0_20),([2],v_0_21),([3],v_0_11),([4],v_0_17)] :: Array [(1,v_0_22)] R) :: P Z
         x_show_0_1 <- discreteUniform 1 v_0_24 :: P Z
         x_show_0_2 <- uniform -0.5 0.5 :: P R
         x_show_0_3 <- discreteUniform 1 input.0_0 :: P Z
         x_show_0_4 <- uniform 0.0 1.0 :: P R
         x_show_0_5 <- uniform 0.0 40907.0 :: P R
         x_show_0_6 <- normal 0.0 1.0 :: P R
         x_show_0_7 <- discreteUniform 1 input.0_0 :: P Z
         x_show_0_8 <- bernoulli v_0_450 :: P B
         return [v_0_452,v_0_453,v_0_454,v_0_455]
In [10]:
samples <- iterateLimit 10 (runStep kernel) (n0,s0,g0, list coalData')
[(n,s,g) | (n,s,g,y) <- samples]
[(4,[5103.668914593937,25700.783608506343,29470.225091303586,36244.890506445045],[4.892522331906208e-3,1.4197347319748954e-2,2.00712577711841e-3,1.3410980131244538e-3,2.615629151408432e-3]),(3,[25700.783608506343,29470.225091303586,36244.890506445045],[4.892522331906208e-3,2.00712577711841e-3,1.3410980131244538e-3,2.615629151408432e-3]),(4,[14899.056729097752,25700.783608506343,29470.225091303586,36244.890506445045],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.3410980131244538e-3,2.615629151408432e-3]),(3,[14899.056729097752,25700.783608506343,29470.225091303586],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.3410980131244538e-3]),(3,[14899.056729097752,25700.783608506343,29470.225091303586],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.3410980131244538e-3]),(3,[14899.056729097752,25700.783608506343,29470.225091303586],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.3410980131244538e-3]),(3,[14899.056729097752,25700.783608506343,29470.225091303586],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.491834168992973e-3]),(3,[14899.056729097752,25700.783608506343,29470.225091303586],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.491834168992973e-3]),(3,[14899.056729097752,25700.783608506343,30337.36331137781],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.491834168992973e-3]),(3,[14899.056729097752,25700.783608506343,30337.36331137781],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3,1.491834168992973e-3]),(2,[14899.056729097752,25700.783608506343],[4.892522331906208e-3,2.413468786837807e-3,2.00712577711841e-3])]
In [11]:
-- sampling takes a long time, so we generate samples offline
import System.IO
output <- readFile "rj-coal.log" -- from stochaskell/app/rj
let samples = read <$> lines output :: [(Int,[Double],[Double])]
In [12]:
take 10 samples
[(3,[24969.58172348803,26526.26524169898,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,4.5984983765520634e-4,2.641155121557546e-3]),(3,[24969.58172348803,26526.26524169898,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,5.601448084736443e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3]),(2,[24969.58172348803,28059.88881142932],[3.94532773540886e-3,3.3475125679093145e-4,2.641155121557546e-3])]
In [13]:
let x = linspace (0, last coalData') 512
    ys = [staircase s g <$> x | (n,s,g) <- samples] :: [[Double]]
    ymean = mean ys
    ymean2 = mean $ square <$> ys
    ystd = sqrt $ ymean2 - square ymean
In [14]:
:opt svg
import Language.Stochaskell.Plot
In [15]:
let x' = (/10000) <$> x
pMain = do
  plot . return $ PlotFillBetween "" (solidFillStyle $ blue `withOpacity` 0.1)
    (zip x' $ zip ((10000 *) <$> ymean + ystd) ((10000 *) <$> ymean - ystd))
  plot $ line "" [zip x' ((10000 *) <$> ymean)]
  xlim (0,4)
  ylim (0,100)
  xlabel "Time (10⁴ days)"
  ylabel "Posterior mean rate"
  layout_margin .= 0
toRenderable pMain
In [16]:
pHist = do
  plot . return . histToPlot $ defaultNormedPlotHist
    { _plot_hist_title = "", _plot_hist_values = [real n | (n,s,g) <- samples], _plot_hist_range = Just (1,8), _plot_hist_bins = 7 }
  xlim (1,8)
  layout_title .= "Number of change points"
  layout_title_style . font_size .= 10
  layout_title_style . font_weight .= FontWeightNormal
  layout_margin .= 0
  layout_top_axis_visibility    .= AxisVisibility True True True
  layout_bottom_axis_visibility .= AxisVisibility True True False
  layout_right_axis_visibility  .= AxisVisibility True True True
  layout_left_axis_visibility   .= AxisVisibility True True False
  layout_x_axis . laxis_generate .= autoScaledAxis def { _la_nLabels = 7, _la_nTicks = 7 }
  layout_y_axis . laxis_generate .= autoScaledAxis def { _la_nLabels = 3, _la_nTicks = 30 }
toRenderable pHist
In [17]:
pX = do
  setColors $ (`withOpacity` 0.75) <$> [blue,green,red,orange,yellow,violet]
  plot $ line "n=1" [kde 625 [s !! (i-1) | (n,s,g) <- samples, n == 1] | i <- [1..1]]
  plot $ line "n=2" [kde 625 [s !! (i-1) | (n,s,g) <- samples, n == 2] | i <- [1..2]]
  plot $ line "n=3" [kde 625 [s !! (i-1) | (n,s,g) <- samples, n == 3] | i <- [1..3]]
  xlim (0,40000)
  layout_legend .= Nothing
  layout_margin .= 0
  layout_bottom_axis_visibility .= AxisVisibility False False True
  layout_left_axis_visibility   .= AxisVisibility False False True
  layout_x_axis . laxis_override .= (\d -> d { _axis_ticks = [], _axis_labels = [] })
  layout_y_axis . laxis_override .= (\d -> d { _axis_ticks = [], _axis_labels = [], _axis_grid = [] })
  ylabel "Change points"
toRenderable pX
In [18]:
import Data.Tuple
pY = do
  setColors $ (`withOpacity` 0.75) <$> [blue,green,red,orange,yellow,violet]
  plot $ line "n=1" [map swap $ kde 0.0003 [g !! i | (n,s,g) <- samples, n == 1] | i <- [0..1]]
  plot $ line "n=2" [map swap $ kde 0.0003 [g !! i | (n,s,g) <- samples, n == 2] | i <- [0..2]]
  plot $ line "n=3" [map swap $ kde 0.0003 [g !! i | (n,s,g) <- samples, n == 3] | i <- [0..3]]
  ylim (0,0.00999)
  layout_legend .= Nothing
  layout_margin .= 0
  layout_bottom_axis_visibility .= AxisVisibility False False True
  layout_left_axis_visibility   .= AxisVisibility False False True
  layout_x_axis . laxis_override .= (\d -> d { _axis_ticks = [], _axis_labels = [], _axis_grid = [] })
  layout_y_axis . laxis_override .= (\d -> d { _axis_ticks = [], _axis_labels = [] })
  xlabel "Step heights"
toRenderable pY
In [19]:
let cell = layoutToGrid . execEC
toRenderable $ (cell pX `beside` cell pHist) `above` (cell pMain `beside` cell pY)