In [1]:
{-# 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
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)))