In [1]:
{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax, DeriveGeneric, DeriveAnyClass #-}

Stochaskell, version 0.1.0
Copyright (C) 2015-2020 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]:
soccer1 :: Z -> P (R,ZVec)
soccer1 n = do
lam <- gamma 25 10
y <- joint vector [ poisson lam
| _ <- 1...n ]
return (lam,y)

soccer2 :: Z -> P (R,R,ZVec)
soccer2 n = do
lam <- gamma 25 10
kap <- gamma 1 10
let a = 1/kap
b = a/lam
y <- joint vector [ negBinomial a b
| _ <- 1...n ]
return (lam,kap,y)

In [3]:
import Language.Stochaskell.Expression
import GHC.Generics (Generic)

In [4]:
data Model = Model1 R ZVec
| Model2 R R ZVec
deriving (Show, Generic, ExprType, Constructor)

In [5]:
soccer :: Z -> P (Expr Model)
soccer n = do
(lam1,     y1) <- soccer1 n
(lam2,kap2,y2) <- soccer2 n
k <- bernoulli 0.5
return $if k then fromConcrete (Model1 lam1 y1) else fromConcrete (Model2 lam2 kap2 y2)  In [6]: soccerJump :: Model -> P (Expr Model) soccerJump (Model1 lam y) = do u <- normal 0 1.5 let kap = 0.015 * exp u return$ fromConcrete (Model2 lam kap y)
soccerJump (Model2 lam _ y) =
return $fromConcrete (Model1 lam y)  In [7]: soccerHMC :: Z -> Model -> IO Model soccerHMC n (Model1 lam y) = do let posterior = [ lam' | (lam',y') <- soccer1 n , y' == y ] samples <- hmcStanInit 10 posterior lam let lam' = last samples return (Model1 lam' y) soccerHMC n (Model2 lam kap y) = do let posterior = [ (lam',kap') | (lam',kap',y') <- soccer2 n , y' == y ] let s0 = (lam,kap) samples <- hmcStanInit 10 posterior s0 let (lam',kap') = last samples return (Model2 lam' kap' y)  In [8]: soccerStep :: Z -> Model -> IO Model soccerStep n m = do m' <- soccerHMC n m m'' <- soccer n rjmcC soccerJump runCC fromConcrete m' return$ fromRight (eval' m'')

In [9]:
-- https://cran.r-project.org/web/packages/engsoccerdata/README.html
totgoal = [5, 1, 3, 2, 3, 2, 2, 5, 3, 1, 0, 7, 2, 4, 4, 2, 4, 5, 6, 0, 4, 1, 4, 1, 2, 4, 0, 2, 1, 2, 5, 3, 1, 3, 2, 0, 3, 2, 2, 1, 3, 1, 1, 0, 1, 1, 4, 3, 4, 3, 0, 5, 1, 2, 2, 3, 2, 1, 2, 2, 0, 5, 1, 2, 3, 1, 2, 7, 5, 3, 3, 2, 0, 2, 5, 2, 2, 2, 1, 0, 5, 2, 1, 3, 4, 2, 3, 2, 2, 1, 2, 1, 2, 5, 2, 1, 0, 2, 2, 1, 2, 0, 2, 2, 7, 4, 3, 4, 3, 2, 5, 0, 2, 1, 1, 3, 2, 6, 6, 2, 3, 5, 2, 2, 3, 1, 3, 2, 2, 3, 4, 5, 1, 1, 5, 0, 1, 4, 4, 2, 4, 4, 1, 2, 1, 1, 1, 4, 1, 4, 3, 1, 4, 6, 0, 3, 3, 3, 1, 1, 2, 3, 5, 1, 1, 4, 3, 1, 7, 3, 1, 1, 4, 2, 1, 1, 0, 5, 4, 6, 1, 0, 2, 2, 3, 1, 1, 1, 2, 3, 4, 4, 5, 0, 1, 5, 1, 2, 3, 1, 4, 1, 3, 3, 3, 2, 0, 3, 1, 2, 1, 3, 3, 5, 4, 1, 2, 6, 1, 2, 0, 2, 3, 0, 2, 3, 1, 4, 3, 4, 1, 2, 7, 3, 3, 1, 5, 0, 0, 5, 3, 2, 2, 6, 4, 2, 5, 1, 2, 1, 1, 4, 0, 1, 2, 2, 4, 1, 2, 4, 2, 5, 4, 3, 0, 4, 2, 2, 2, 4, 2, 3, 2, 1, 1, 4, 3, 4, 1, 0, 3, 2, 1, 2, 2, 3, 4, 1, 1, 0, 4, 3, 1, 3, 2, 3, 4, 3, 5, 5, 2, 2, 2, 1, 2, 0, 2, 5, 1, 4, 2, 2, 1, 0, 3, 3, 2, 2, 4, 5, 3, 2, 4, 3, 3, 5, 2, 0, 3, 3, 4, 0, 2, 2, 3, 2, 3, 3, 1, 2, 1, 3, 0, 4, 3, 4, 3, 0, 4, 4, 3, 3, 1, 3, 3, 6, 6, 2, 3, 1, 2, 5, 5, 2, 3, 3, 3, 1, 2, 1, 1, 7, 3, 2, 1, 3, 1, 3, 1, 3, 2, 8, 3, 4, 2, 2, 4, 3, 4, 3, 2, 2, 4, 3, 3, 3, 3, 1, 3, 1, 2, 1, 2, 0, 2, 2, 0, 4, 3, 2, 2, 0, 3, 3, 2, 2, 1, 2, 2, 3, 1, 5, 2, 2, 2, 1, 6, 1, 3, 4, 3, 6, 3, 2, 4, 3, 3, 4, 4, 3, 2, 1, 2, 3, 2, 0, 4, 0, 3, 5, 4, 1, 2, 1, 4, 1, 3, 3, 1, 2, 1, 2, 4, 3, 1, 3, 4, 2, 1, 0, 2, 2, 0, 4, 1, 2, 2, 3, 4, 3, 2, 4, 1, 3, 0, 3, 1, 3, 4, 3, 1, 4, 1, 4, 1, 1, 1, 1, 3, 5, 5, 3, 2, 6, 0, 3, 3, 2, 2, 3, 3, 2, 4, 3, 2, 2, 2, 3, 2, 1, 1, 4, 3, 3, 3, 2, 1, 1, 2, 0, 0, 1, 5, 4, 2, 3, 4, 2, 0, 4, 1, 1, 2, 2, 0, 2, 4, 3, 2, 3, 2, 1, 2, 3, 2, 0, 1, 3, 4, 0, 1, 1, 0, 0, 2, 0, 3, 0, 2, 1, 1, 4, 5, 5, 2, 2, 3, 6, 2, 4, 2, 2, 3, 5, 2, 1, 4, 1, 4, 2, 4, 1, 6, 2, 3, 3, 4, 0, 2, 3, 1, 4, 3, 4, 5, 5, 1, 2, 0, 4, 2, 3, 0, 0, 2, 3, 3, 1, 4, 0, 1, 5, 1, 4, 3, 4, 3, 0, 4, 3, 1, 1, 2, 2, 2, 3, 3, 3, 0, 3, 4, 4, 2, 3, 2, 1, 4, 2, 3, 1, 2, 1, 2, 1, 3, 1, 2, 5, 1, 0, 4, 4, 2, 6, 5, 1, 4, 0, 4, 3, 2, 2, 2, 2, 1, 3, 3, 3, 2, 3, 3, 1, 3, 3, 4, 3, 2, 5, 6, 3, 2, 0, 1, 3, 4, 3, 5, 3, 1, 2, 4, 1, 4, 3, 0, 3, 1, 4, 1, 3, 6, 3, 2, 3, 2, 2, 6, 0, 1, 0, 2, 2, 1, 2, 3, 4, 4, 5, 1, 6, 3, 1, 1, 2, 2, 3, 1, 1, 7, 1, 2, 1, 0, 3, 4, 5, 5, 2, 0, 4, 4, 4, 1, 1, 1, 1, 1, 6, 2, 3, 2, 2, 2, 2, 1, 5, 1, 3, 2, 1, 4, 2, 3, 4, 2, 5, 3, 2, 2, 3, 6, 2, 4, 2, 2, 2, 3, 3, 2, 5, 2, 5, 4, 4, 1, 3, 1, 2, 4, 3, 5, 1, 1, 2, 2, 2, 4, 4, 1, 3, 2, 2, 2, 4, 5, 1, 5, 2, 4, 3, 5, 1, 4, 0, 2, 0, 1, 2, 2, 4, 1, 6, 1, 2, 1, 4, 5, 2, 3, 3, 1, 1, 3, 0, 4, 0, 1, 0, 4, 1, 3, 2, 2, 1, 5, 3, 8, 5, 0, 2, 7, 2, 0, 0, 6, 3, 1, 3, 1, 1, 2, 2, 1, 2, 8, 6, 3, 3, 2, 2, 2, 4, 3, 2, 1, 1, 1, 4, 4, 0, 3, 5, 1, 5, 4, 4, 2, 2, 1, 1, 3, 3, 1, 1, 2, 4, 4, 1, 8, 0, 2, 3, 3, 3, 2, 4, 3, 3, 0, 1, 2, 6, 3, 3, 1, 2, 4, 4, 6, 1, 2, 2, 4, 0, 4, 4, 2, 6, 1, 2, 1, 1, 5, 3, 5, 3, 3, 4, 4, 2, 4, 1, 1, 4, 6, 2, 1, 2, 5, 0, 1, 4, 4, 4, 3, 1, 3, 2, 0, 3, 4, 1, 2, 2, 2, 5, 3, 2, 3, 3, 5, 6, 2, 0, 1, 1, 5, 4, 3, 3, 2, 3, 1, 2, 1, 2, 1, 2, 9, 4, 4, 2, 1, 4, 2, 3, 1, 2, 0, 3, 1, 0, 2, 4, 5, 2, 3, 2, 6, 2, 5, 3, 2, 4, 4, 1, 0, 2, 6, 1, 4, 2, 4, 0, 1, 0, 0, 2, 1, 0, 11, 1, 1, 0, 2, 4, 3, 3, 0, 2, 3, 1, 1, 2, 4, 2, 2, 2, 3, 2, 3, 1, 3, 3, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 5, 2, 2, 3, 1, 3, 2, 4, 8, 5, 3, 2, 8, 4, 4, 6, 2, 3, 2, 2, 5, 2, 10, 2, 4, 4, 1, 4, 2, 3, 2, 4, 3, 2, 3, 1, 2, 3, 3, 4, 1, 2, 4, 2, 2, 0, 3, 2, 8, 1, 2, 2, 3, 2, 1, 2, 2, 1, 1, 2, 0, 3, 2, 1]

In [10]:
let goalData = list totgoal
n = integer (length totgoal)

In [11]:
compileCC (soccer n rjmcC soccerJump)

/*
let v_0_0 = getExternal x_cc_0_0 :: Union[(R, ZVec), (R, R, ZVec)]
in do x_cc_0_0 <- switch i_0_0 of
C0 i_1_0 i_1_1 ->
let v_1_0 = getExternal x_cc_1_0 :: R
v_1_1 = gamma_lpdf i_1_0 25.0 10.0 :: R
v_1_2 = exp x_cc_1_0 :: R
v_1_3 = 1.5e-2 * v_1_2 :: R
v_1_4 = gamma_lpdf v_1_3 1.0 10.0 :: R
v_1_5 = 1.0 / v_1_3 :: R
v_1_6 = v_1_5 / i_1_0 :: R
v_1_7 = vectorSize i_1_1 :: Z
v_1_8 =
[ i_2_1
| i_2_1 <- 1...v_1_7 ] :: ZVec
v_1_9 = foldl 0.0 v_1_8 $\i_2_12 i_2_11 -> let v_2_0 = neg_binomial_lpdf i_1_1!i_2_11 v_1_5 v_1_6 :: R v_2_1 = i_2_12 + v_2_0 :: R in v_2_1 :: R v_1_10 = +s -0.6931471805599453 v_1_1 v_1_4 v_1_9 :: R v_1_11 = foldl 0.0 v_1_8$ \i_2_12 i_2_11 ->
let v_2_0 = poisson_lpdf i_1_1!i_2_11 i_1_0 :: R
v_2_1 = i_2_12 + v_2_0 :: R
in v_2_1 :: R
v_1_12 = +s -0.6931471805599453 v_1_1 v_1_11 :: R
v_1_13 = v_1_10 - v_1_12 :: R
v_1_14 = normal_lpdf x_cc_1_0 0.0 1.5 :: R
v_1_15 = negate v_1_14 :: R
v_1_16 = +s -4.199705077879927 x_cc_1_0 v_1_15 :: R
v_1_17 = +s x_cc_1_0 v_1_13 v_1_15 -4.199705077879927 :: R
v_1_18 = exp v_1_17 :: R
v_1_19 = min 1.0 v_1_18 :: R
v_1_20 = getExternal x_cc_1_1 :: B
v_1_21 = ifThenElse x_cc_1_1 (C1 i_1_0 v_1_3 i_1_1 :: Union[(R, ZVec), (R, R, ZVec)]) (C0 i_1_0 i_1_1 :: Union[(R, ZVec), (R, R, ZVec)]) :: Union[(R, ZVec), (R, R, ZVec)]
in do x_cc_1_0 <- normal 0.0 1.5 :: P R
x_cc_1_1 <- bernoulli v_1_19 :: P B
return [v_1_21]
C1 i_1_0 i_1_1 i_1_2 ->
let v_1_0 = gamma_lpdf i_1_0 25.0 10.0 :: R
v_1_1 = vectorSize i_1_2 :: Z
v_1_2 =
[ i_2_1
| i_2_1 <- 1...v_1_1 ] :: ZVec
v_1_3 = foldl 0.0 v_1_2 $\i_2_12 i_2_11 -> let v_2_0 = poisson_lpdf i_1_2!i_2_11 i_1_0 :: R v_2_1 = i_2_12 + v_2_0 :: R in v_2_1 :: R v_1_4 = +s -0.6931471805599453 v_1_0 v_1_3 :: R v_1_5 = gamma_lpdf i_1_1 1.0 10.0 :: R v_1_6 = 1.0 / i_1_1 :: R v_1_7 = v_1_6 / i_1_0 :: R v_1_8 = foldl 0.0 v_1_2$ \i_2_12 i_2_11 ->
let v_2_0 = neg_binomial_lpdf i_1_2!i_2_11 v_1_6 v_1_7 :: R
v_2_1 = i_2_12 + v_2_0 :: R
in v_2_1 :: R
v_1_9 = +s -0.6931471805599453 v_1_0 v_1_5 v_1_8 :: R
v_1_10 = v_1_4 - v_1_9 :: R
v_1_11 = log i_1_1 :: R
v_1_12 = v_1_11 - -4.199705077879927 :: R
v_1_13 = normal_lpdf v_1_12 0.0 1.5 :: R
v_1_14 = negate v_1_11 :: R
v_1_15 = v_1_13 + v_1_14 :: R
v_1_16 = +s v_1_10 v_1_13 v_1_14 :: R
v_1_17 = exp v_1_16 :: R
v_1_18 = min 1.0 v_1_17 :: R
v_1_19 = getExternal x_cc_1_0 :: B
v_1_20 = ifThenElse x_cc_1_0 (C0 i_1_0 i_1_2 :: Union[(R, ZVec), (R, R, ZVec)]) (C1 i_1_0 i_1_1 i_1_2 :: Union[(R, ZVec), (R, R, ZVec)]) :: Union[(R, ZVec), (R, R, ZVec)]
in do x_cc_1_0 <- bernoulli v_1_18 :: P B
return [v_1_20]
return [x_cc_0_0]
*/

#include <cmath>
#include <iostream>
#include <random>
#include <tuple>
#include <boost/math/distributions.hpp>
#include <boost/math/special_functions/factorials.hpp>
#define eigen_assert(X) do { if(!(X)) throw std::runtime_error(#X); } while(false);
#include <Eigen/Dense>
#include <mpark/variant.hpp>
#include <backward.hpp>
using namespace std;
using namespace Eigen;
using namespace mpark;
namespace backward { backward::SignalHandling sh; }
IOFormat eigenFormat(StreamPrecision, DontAlignCols, ",", ";", "", "", "[", "]");
#define LOG_DET(m) ((m).isLowerTriangular() ? (m).diagonal().array().log().sum() : (m).llt().matrixL().toDenseMatrix().diagonal().array().log().sum() * 2)

int main() {
random_device rd;
mt19937 gen(rd());
variant<tuple<double, Matrix<int, Dynamic, 1>>, tuple<double, double, Matrix<int, Dynamic, 1>>> i_0_0; int i_1_0; cin >> i_1_0;
switch(i_1_0) {
case 0: {
double i_1_1; cin >> i_1_1;
Matrix<int, Dynamic, 1> i_1_2; int i_1_2_length; cin >> i_1_2_length;
i_1_2.resize(i_1_2_length);
for(int i_2_1 = 1; i_2_1 <= i_1_2_length; i_2_1++) {
cin >> i_1_2(i_2_1-1);
}
i_0_0 = make_tuple(i_1_1, i_1_2);
break;
}
case 1: {
double i_1_1; cin >> i_1_1;
double i_1_2; cin >> i_1_2;
Matrix<int, Dynamic, 1> i_1_3; int i_1_3_length; cin >> i_1_3_length;
i_1_3.resize(i_1_3_length);
for(int i_2_1 = 1; i_2_1 <= i_1_3_length; i_2_1++) {
cin >> i_1_3(i_2_1-1);
}
i_0_0 = make_tuple(i_1_1, i_1_2, i_1_3);
break;
}
}
variant<tuple<double, Matrix<int, Dynamic, 1>>, tuple<double, double, Matrix<int, Dynamic, 1>>> x_cc_0_0; try { switch(i_0_0.index()) {
case 0: {
double i_1_0; Matrix<int, Dynamic, 1> i_1_1; tie(i_1_0, i_1_1) = get<0>(i_0_0);
double x_cc_1_0; try { x_cc_1_0 = std::normal_distribution<>(0, 1.5)(gen); } catch(const std::runtime_error& e) { cerr << "x_cc_1_0: " << e.what() << endl; }
double v_1_1; try { v_1_1 = log(boost::math::pdf(boost::math::gamma_distribution<>(25, 1./10), i_1_0)); } catch(const std::runtime_error& e) { cerr << "v_1_1: " << e.what() << endl; }
double v_1_2; try { v_1_2 = exp(x_cc_1_0); } catch(const std::runtime_error& e) { cerr << "v_1_2: " << e.what() << endl; }
double v_1_3; try { v_1_3 = 1.5e-2 * v_1_2; } catch(const std::runtime_error& e) { cerr << "v_1_3: " << e.what() << endl; }
double v_1_4; try { v_1_4 = log(boost::math::pdf(boost::math::gamma_distribution<>(1, 1./10), v_1_3)); } catch(const std::runtime_error& e) { cerr << "v_1_4: " << e.what() << endl; }
double v_1_5; try { v_1_5 = (double) 1 / (double) v_1_3; } catch(const std::runtime_error& e) { cerr << "v_1_5: " << e.what() << endl; }
double v_1_6; try { v_1_6 = (double) v_1_5 / (double) i_1_0; } catch(const std::runtime_error& e) { cerr << "v_1_6: " << e.what() << endl; }
int v_1_7; try { v_1_7 = i_1_1.size(); } catch(const std::runtime_error& e) { cerr << "v_1_7: " << e.what() << endl; }
Matrix<int, Dynamic, 1> v_1_8; try { v_1_8.resize(v_1_7);
for(int i_2_1 = 1; i_2_1 <= v_1_7; i_2_1++) {

v_1_8(i_2_1-1) = i_2_1;
} } catch(const std::runtime_error& e) { cerr << "v_1_8: " << e.what() << endl; }
double v_1_9; try { v_1_9 = 0;
for(int i_2_0 = 1; i_2_0 <= v_1_8.size(); i_2_0++) {
int i_2_11 = v_1_8[i_2_0-1];
double i_2_12 = v_1_9;
double v_2_0; v_2_0 = log(boost::math::pdf(boost::math::negative_binomial_distribution<>(v_1_5, v_1_6 / (v_1_6 + 1)), i_1_1(i_2_11-1)));
double v_2_1; v_2_1 = i_2_12 + v_2_0;
v_1_9 = v_2_1;

} } catch(const std::runtime_error& e) { cerr << "v_1_9: " << e.what() << endl; }
double v_1_10; try { v_1_10 = -0.6931471805599453 + v_1_1 + v_1_4 + v_1_9; } catch(const std::runtime_error& e) { cerr << "v_1_10: " << e.what() << endl; }
double v_1_11; try { v_1_11 = 0;
for(int i_2_0 = 1; i_2_0 <= v_1_8.size(); i_2_0++) {
int i_2_11 = v_1_8[i_2_0-1];
double i_2_12 = v_1_11;
double v_2_0; v_2_0 = log(boost::math::pdf(boost::math::poisson_distribution<>(i_1_0), i_1_1(i_2_11-1)));
double v_2_1; v_2_1 = i_2_12 + v_2_0;
v_1_11 = v_2_1;

} } catch(const std::runtime_error& e) { cerr << "v_1_11: " << e.what() << endl; }
double v_1_12; try { v_1_12 = -0.6931471805599453 + v_1_1 + v_1_11; } catch(const std::runtime_error& e) { cerr << "v_1_12: " << e.what() << endl; }
double v_1_13; try { v_1_13 = v_1_10 - v_1_12; } catch(const std::runtime_error& e) { cerr << "v_1_13: " << e.what() << endl; }
double v_1_14; try { v_1_14 = log(boost::math::pdf(boost::math::normal_distribution<>(0, 1.5), x_cc_1_0)); } catch(const std::runtime_error& e) { cerr << "v_1_14: " << e.what() << endl; }
double v_1_15; try { v_1_15 = -v_1_14; } catch(const std::runtime_error& e) { cerr << "v_1_15: " << e.what() << endl; }
double v_1_16; try { v_1_16 = -4.199705077879927 + x_cc_1_0 + v_1_15; } catch(const std::runtime_error& e) { cerr << "v_1_16: " << e.what() << endl; }
double v_1_17; try { v_1_17 = x_cc_1_0 + v_1_13 + v_1_15 + -4.199705077879927; } catch(const std::runtime_error& e) { cerr << "v_1_17: " << e.what() << endl; }
double v_1_18; try { v_1_18 = exp(v_1_17); } catch(const std::runtime_error& e) { cerr << "v_1_18: " << e.what() << endl; }
double v_1_19; try { v_1_19 = min<double>(1, v_1_18); } catch(const std::runtime_error& e) { cerr << "v_1_19: " << e.what() << endl; }
int x_cc_1_1; try { do { x_cc_1_1 = std::bernoulli_distribution(v_1_19)(gen); } while(x_cc_1_1 < 0 || x_cc_1_1 > 1); } catch(const std::runtime_error& e) { cerr << "x_cc_1_1: " << e.what() << endl; }
variant<tuple<double, Matrix<int, Dynamic, 1>>, tuple<double, double, Matrix<int, Dynamic, 1>>> v_1_21; try { if(x_cc_1_1) v_1_21 = make_tuple(i_1_0, v_1_3, i_1_1); else v_1_21 = make_tuple(i_1_0, i_1_1); } catch(const std::runtime_error& e) { cerr << "v_1_21: " << e.what() << endl; }
x_cc_0_0 = v_1_21; break;
}
case 1: {
double i_1_0; double i_1_1; Matrix<int, Dynamic, 1> i_1_2; tie(i_1_0, i_1_1, i_1_2) = get<1>(i_0_0);
double v_1_0; try { v_1_0 = log(boost::math::pdf(boost::math::gamma_distribution<>(25, 1./10), i_1_0)); } catch(const std::runtime_error& e) { cerr << "v_1_0: " << e.what() << endl; }
int v_1_1; try { v_1_1 = i_1_2.size(); } catch(const std::runtime_error& e) { cerr << "v_1_1: " << e.what() << endl; }
Matrix<int, Dynamic, 1> v_1_2; try { v_1_2.resize(v_1_1);
for(int i_2_1 = 1; i_2_1 <= v_1_1; i_2_1++) {

v_1_2(i_2_1-1) = i_2_1;
} } catch(const std::runtime_error& e) { cerr << "v_1_2: " << e.what() << endl; }
double v_1_3; try { v_1_3 = 0;
for(int i_2_0 = 1; i_2_0 <= v_1_2.size(); i_2_0++) {
int i_2_11 = v_1_2[i_2_0-1];
double i_2_12 = v_1_3;
double v_2_0; v_2_0 = log(boost::math::pdf(boost::math::poisson_distribution<>(i_1_0), i_1_2(i_2_11-1)));
double v_2_1; v_2_1 = i_2_12 + v_2_0;
v_1_3 = v_2_1;

} } catch(const std::runtime_error& e) { cerr << "v_1_3: " << e.what() << endl; }
double v_1_4; try { v_1_4 = -0.6931471805599453 + v_1_0 + v_1_3; } catch(const std::runtime_error& e) { cerr << "v_1_4: " << e.what() << endl; }
double v_1_5; try { v_1_5 = log(boost::math::pdf(boost::math::gamma_distribution<>(1, 1./10), i_1_1)); } catch(const std::runtime_error& e) { cerr << "v_1_5: " << e.what() << endl; }
double v_1_6; try { v_1_6 = (double) 1 / (double) i_1_1; } catch(const std::runtime_error& e) { cerr << "v_1_6: " << e.what() << endl; }
double v_1_7; try { v_1_7 = (double) v_1_6 / (double) i_1_0; } catch(const std::runtime_error& e) { cerr << "v_1_7: " << e.what() << endl; }
double v_1_8; try { v_1_8 = 0;
for(int i_2_0 = 1; i_2_0 <= v_1_2.size(); i_2_0++) {
int i_2_11 = v_1_2[i_2_0-1];
double i_2_12 = v_1_8;
double v_2_0; v_2_0 = log(boost::math::pdf(boost::math::negative_binomial_distribution<>(v_1_6, v_1_7 / (v_1_7 + 1)), i_1_2(i_2_11-1)));
double v_2_1; v_2_1 = i_2_12 + v_2_0;
v_1_8 = v_2_1;

} } catch(const std::runtime_error& e) { cerr << "v_1_8: " << e.what() << endl; }
double v_1_9; try { v_1_9 = -0.6931471805599453 + v_1_0 + v_1_5 + v_1_8; } catch(const std::runtime_error& e) { cerr << "v_1_9: " << e.what() << endl; }
double v_1_10; try { v_1_10 = v_1_4 - v_1_9; } catch(const std::runtime_error& e) { cerr << "v_1_10: " << e.what() << endl; }
double v_1_11; try { v_1_11 = log(i_1_1); } catch(const std::runtime_error& e) { cerr << "v_1_11: " << e.what() << endl; }
double v_1_12; try { v_1_12 = v_1_11 - -4.199705077879927; } catch(const std::runtime_error& e) { cerr << "v_1_12: " << e.what() << endl; }
double v_1_13; try { v_1_13 = log(boost::math::pdf(boost::math::normal_distribution<>(0, 1.5), v_1_12)); } catch(const std::runtime_error& e) { cerr << "v_1_13: " << e.what() << endl; }
double v_1_14; try { v_1_14 = -v_1_11; } catch(const std::runtime_error& e) { cerr << "v_1_14: " << e.what() << endl; }
double v_1_15; try { v_1_15 = v_1_13 + v_1_14; } catch(const std::runtime_error& e) { cerr << "v_1_15: " << e.what() << endl; }
double v_1_16; try { v_1_16 = v_1_10 + v_1_13 + v_1_14; } catch(const std::runtime_error& e) { cerr << "v_1_16: " << e.what() << endl; }
double v_1_17; try { v_1_17 = exp(v_1_16); } catch(const std::runtime_error& e) { cerr << "v_1_17: " << e.what() << endl; }
double v_1_18; try { v_1_18 = min<double>(1, v_1_17); } catch(const std::runtime_error& e) { cerr << "v_1_18: " << e.what() << endl; }
int x_cc_1_0; try { do { x_cc_1_0 = std::bernoulli_distribution(v_1_18)(gen); } while(x_cc_1_0 < 0 || x_cc_1_0 > 1); } catch(const std::runtime_error& e) { cerr << "x_cc_1_0: " << e.what() << endl; }
variant<tuple<double, Matrix<int, Dynamic, 1>>, tuple<double, double, Matrix<int, Dynamic, 1>>> v_1_20; try { if(x_cc_1_0) v_1_20 = make_tuple(i_1_0, i_1_2); else v_1_20 = make_tuple(i_1_0, i_1_1, i_1_2); } catch(const std::runtime_error& e) { cerr << "v_1_20: " << e.what() << endl; }
x_cc_0_0 = v_1_20; break;
}
} } catch(const std::runtime_error& e) { cerr << "x_cc_0_0: " << e.what() << endl; }
switch(x_cc_0_0.index()) {
case 0: {
double i_1_1; Matrix<int, Dynamic, 1> i_1_2; tie(i_1_1, i_1_2) = get<0>(x_cc_0_0);
cout << 0 << ' ';
cout << i_1_1 << ' ';
cout << i_1_2.size() << ' ';
for(int i_2_1 = 1; i_2_1 <= i_1_2.size(); i_2_1++) {
cout << i_1_2(i_2_1-1) << ' ';
}
break;
}
case 1: {
double i_1_1; double i_1_2; Matrix<int, Dynamic, 1> i_1_3; tie(i_1_1, i_1_2, i_1_3) = get<1>(x_cc_0_0);
cout << 1 << ' ';
cout << i_1_1 << ' ';
cout << i_1_2 << ' ';
cout << i_1_3.size() << ' ';
for(int i_2_1 = 1; i_2_1 <= i_1_3.size(); i_2_1++) {
cout << i_1_3(i_2_1-1) << ' ';
}
break;
}
} cout << endl;
}

g++ -std=c++11 -Ibackward-cpp -Ieigen -Ivariant/include -DBACKWARD_HAS_BFD=1 -g -rdynamic /home/jovyan/stochaskell/cache/cc/sampler_576e8212a8171ab5532cf51f6363c4521cf709eb.cc -lbfd -ldl -o /home/jovyan/stochaskell/cache/cc/sampler_576e8212a8171ab5532cf51f6363c4521cf709eb
In [12]:
-- reduce number of iterations from 5000 to speed this up
samples <- silence' $iterateLimit 5000 (soccerStep n) (Model1 1 goalData)  In [13]: putStrLn "Posterior probability in favour of model 1:" sum [case m of Model1{} -> 1; _ -> 0 | m <- samples] / genericLength samples  Posterior probability in favour of model 1: 0.7056588682263547 In [14]: :opt svg import Language.Stochaskell.Plot  In [15]: let plot' = layoutToGrid . execEC . plot hist s vals rng = return . histToPlot$ defaultNormedPlotHist
{ _plot_hist_title = s, _plot_hist_values = vals, _plot_hist_range = Just rng }

toRenderable $(plot' (hist "λ₁" [real lam | Model1 lam _ <- samples] (2.4,2.7)) beside plot' (line "lpdf" [[(i, real$ soccer n lpdfAuxC fromConcrete m) | (i,m) <- [1..] zip tail samples]]))
above
(plot' (hist "λ₂" [real lam | Model2 lam kap _ <- samples] (2.4,2.7)) beside
plot' (hist "κ₂" [real kap | Model2 lam kap _ <- samples] (0,0.08)))