{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax, DeriveGeneric, DeriveAnyClass #-}
import Language.Stochaskell
stochaskell
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
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)
import Language.Stochaskell.Expression
import GHC.Generics (Generic)
data Model = Model1 R ZVec
| Model2 R R ZVec
deriving (Show, Generic, ExprType, Constructor)
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)
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)
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)
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'')
-- 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]
let goalData = list totgoal
n = integer (length totgoal)
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
-- reduce number of iterations from 5000 to speed this up
samples <- silence' $ iterateLimit 5000 (soccerStep n) (Model1 1 goalData)
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
:opt svg
import Language.Stochaskell.Plot
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)))