{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Stochaskell, version 0.1.0\n", "Copyright (C) 2015-2020 David A Roberts\n", "This program comes with ABSOLUTELY NO WARRANTY.\n", "This is free software, and you are welcome to redistribute it\n", "under certain conditions; see the LICENSE for details.\n", "\n", "Using installation directory at \n", " /home/jovyan/stochaskell" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax, DeriveGeneric, DeriveAnyClass #-}\n", "import Language.Stochaskell\n", "stochaskell" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "soccer1 :: Z -> P (R,ZVec)\n", "soccer1 n = do\n", " lam <- gamma 25 10\n", " y <- joint vector [ poisson lam\n", " | _ <- 1...n ]\n", " return (lam,y)\n", "\n", "soccer2 :: Z -> P (R,R,ZVec)\n", "soccer2 n = do\n", " lam <- gamma 25 10\n", " kap <- gamma 1 10\n", " let a = 1/kap\n", " b = a/lam\n", " y <- joint vector [ negBinomial a b\n", " | _ <- 1...n ]\n", " return (lam,kap,y)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import Language.Stochaskell.Expression\n", "import GHC.Generics (Generic)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "data Model = Model1 R ZVec\n", " | Model2 R R ZVec\n", " deriving (Show, Generic, ExprType, Constructor)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "soccer :: Z -> P (Expr Model)\n", "soccer n = do\n", " (lam1, y1) <- soccer1 n\n", " (lam2,kap2,y2) <- soccer2 n\n", " k <- bernoulli 0.5\n", " return $ if k\n", " then fromConcrete (Model1 lam1 y1)\n", " else fromConcrete (Model2 lam2 kap2 y2)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "soccerJump :: Model -> P (Expr Model)\n", "soccerJump (Model1 lam y) = do\n", " u <- normal 0 1.5\n", " let kap = 0.015 * exp u\n", " return $ fromConcrete (Model2 lam kap y)\n", "soccerJump (Model2 lam _ y) =\n", " return $ fromConcrete (Model1 lam y)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "soccerHMC :: Z -> Model -> IO Model\n", "soccerHMC n (Model1 lam y) = do\n", " let posterior =\n", " [ lam'\n", " | (lam',y') <- soccer1 n\n", " , y' == y ]\n", " samples <- hmcStanInit 10 posterior lam\n", " let lam' = last samples\n", " return (Model1 lam' y)\n", "soccerHMC n (Model2 lam kap y) = do\n", " let posterior =\n", " [ (lam',kap')\n", " | (lam',kap',y') <- soccer2 n\n", " , y' == y ]\n", " let s0 = (lam,kap)\n", " samples <- hmcStanInit 10 posterior s0\n", " let (lam',kap') = last samples\n", " return (Model2 lam' kap' y)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "soccerStep :: Z -> Model -> IO Model\n", "soccerStep n m = do\n", " m' <- soccerHMC n m\n", " m'' <- soccer n `rjmcC` soccerJump\n", " `runCC` fromConcrete m'\n", " return $ fromRight (eval' m'')" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "-- https://cran.r-project.org/web/packages/engsoccerdata/README.html\n", "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]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "let goalData = list totgoal\n", " n = integer (length totgoal)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "/*\n", "let v_0_0 = getExternal x_cc_0_0 :: Union[(R, ZVec), (R, R, ZVec)]\n", " in do x_cc_0_0 <- switch i_0_0 of\n", " C0 i_1_0 i_1_1 ->\n", " let v_1_0 = getExternal x_cc_1_0 :: R\n", " v_1_1 = gamma_lpdf i_1_0 25.0 10.0 :: R\n", " v_1_2 = exp x_cc_1_0 :: R\n", " v_1_3 = 1.5e-2 * v_1_2 :: R\n", " v_1_4 = gamma_lpdf v_1_3 1.0 10.0 :: R\n", " v_1_5 = 1.0 / v_1_3 :: R\n", " v_1_6 = v_1_5 / i_1_0 :: R\n", " v_1_7 = vectorSize i_1_1 :: Z\n", " v_1_8 = \n", " [ i_2_1\n", " | i_2_1 <- 1...v_1_7 ] :: ZVec\n", " v_1_9 = foldl 0.0 v_1_8 $ \\i_2_12 i_2_11 ->\n", " let v_2_0 = neg_binomial_lpdf i_1_1!i_2_11 v_1_5 v_1_6 :: R\n", " v_2_1 = i_2_12 + v_2_0 :: R\n", " in v_2_1 :: R\n", " v_1_10 = +s -0.6931471805599453 v_1_1 v_1_4 v_1_9 :: R\n", " v_1_11 = foldl 0.0 v_1_8 $ \\i_2_12 i_2_11 ->\n", " let v_2_0 = poisson_lpdf i_1_1!i_2_11 i_1_0 :: R\n", " v_2_1 = i_2_12 + v_2_0 :: R\n", " in v_2_1 :: R\n", " v_1_12 = +s -0.6931471805599453 v_1_1 v_1_11 :: R\n", " v_1_13 = v_1_10 - v_1_12 :: R\n", " v_1_14 = normal_lpdf x_cc_1_0 0.0 1.5 :: R\n", " v_1_15 = negate v_1_14 :: R\n", " v_1_16 = +s -4.199705077879927 x_cc_1_0 v_1_15 :: R\n", " v_1_17 = +s x_cc_1_0 v_1_13 v_1_15 -4.199705077879927 :: R\n", " v_1_18 = exp v_1_17 :: R\n", " v_1_19 = min 1.0 v_1_18 :: R\n", " v_1_20 = getExternal x_cc_1_1 :: B\n", " 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)]\n", " in do x_cc_1_0 <- normal 0.0 1.5 :: P R\n", " x_cc_1_1 <- bernoulli v_1_19 :: P B\n", " return [v_1_21]\n", " C1 i_1_0 i_1_1 i_1_2 ->\n", " let v_1_0 = gamma_lpdf i_1_0 25.0 10.0 :: R\n", " v_1_1 = vectorSize i_1_2 :: Z\n", " v_1_2 = \n", " [ i_2_1\n", " | i_2_1 <- 1...v_1_1 ] :: ZVec\n", " v_1_3 = foldl 0.0 v_1_2 $ \\i_2_12 i_2_11 ->\n", " let v_2_0 = poisson_lpdf i_1_2!i_2_11 i_1_0 :: R\n", " v_2_1 = i_2_12 + v_2_0 :: R\n", " in v_2_1 :: R\n", " v_1_4 = +s -0.6931471805599453 v_1_0 v_1_3 :: R\n", " v_1_5 = gamma_lpdf i_1_1 1.0 10.0 :: R\n", " v_1_6 = 1.0 / i_1_1 :: R\n", " v_1_7 = v_1_6 / i_1_0 :: R\n", " v_1_8 = foldl 0.0 v_1_2 $ \\i_2_12 i_2_11 ->\n", " let v_2_0 = neg_binomial_lpdf i_1_2!i_2_11 v_1_6 v_1_7 :: R\n", " v_2_1 = i_2_12 + v_2_0 :: R\n", " in v_2_1 :: R\n", " v_1_9 = +s -0.6931471805599453 v_1_0 v_1_5 v_1_8 :: R\n", " v_1_10 = v_1_4 - v_1_9 :: R\n", " v_1_11 = log i_1_1 :: R\n", " v_1_12 = v_1_11 - -4.199705077879927 :: R\n", " v_1_13 = normal_lpdf v_1_12 0.0 1.5 :: R\n", " v_1_14 = negate v_1_11 :: R\n", " v_1_15 = v_1_13 + v_1_14 :: R\n", " v_1_16 = +s v_1_10 v_1_13 v_1_14 :: R\n", " v_1_17 = exp v_1_16 :: R\n", " v_1_18 = min 1.0 v_1_17 :: R\n", " v_1_19 = getExternal x_cc_1_0 :: B\n", " 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)]\n", " in do x_cc_1_0 <- bernoulli v_1_18 :: P B\n", " return [v_1_20]\n", " return [x_cc_0_0]\n", "*/\n", "\n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#define eigen_assert(X) do { if(!(X)) throw std::runtime_error(#X); } while(false);\n", "#include \n", "#include \n", "#include \n", "using namespace std;\n", "using namespace Eigen;\n", "using namespace mpark;\n", "namespace backward { backward::SignalHandling sh; }\n", "IOFormat eigenFormat(StreamPrecision, DontAlignCols, \",\", \";\", \"\", \"\", \"[\", \"]\");\n", "#define LOG_DET(m) ((m).isLowerTriangular() ? (m).diagonal().array().log().sum() : (m).llt().matrixL().toDenseMatrix().diagonal().array().log().sum() * 2)\n", "\n", "int main() {\n", " random_device rd;\n", " mt19937 gen(rd());\n", " variant>, tuple>> i_0_0; int i_1_0; cin >> i_1_0;\n", " switch(i_1_0) {\n", " case 0: {\n", " double i_1_1; cin >> i_1_1;\n", " Matrix i_1_2; int i_1_2_length; cin >> i_1_2_length;\n", " i_1_2.resize(i_1_2_length);\n", " for(int i_2_1 = 1; i_2_1 <= i_1_2_length; i_2_1++) {\n", " cin >> i_1_2(i_2_1-1);\n", " }\n", " i_0_0 = make_tuple(i_1_1, i_1_2);\n", " break;\n", " }\n", " case 1: {\n", " double i_1_1; cin >> i_1_1;\n", " double i_1_2; cin >> i_1_2;\n", " Matrix i_1_3; int i_1_3_length; cin >> i_1_3_length;\n", " i_1_3.resize(i_1_3_length);\n", " for(int i_2_1 = 1; i_2_1 <= i_1_3_length; i_2_1++) {\n", " cin >> i_1_3(i_2_1-1);\n", " }\n", " i_0_0 = make_tuple(i_1_1, i_1_2, i_1_3);\n", " break;\n", " }\n", " }\n", " variant>, tuple>> x_cc_0_0; try { switch(i_0_0.index()) {\n", " case 0: {\n", " double i_1_0; Matrix i_1_1; tie(i_1_0, i_1_1) = get<0>(i_0_0);\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " Matrix v_1_8; try { v_1_8.resize(v_1_7);\n", " for(int i_2_1 = 1; i_2_1 <= v_1_7; i_2_1++) {\n", " \n", " v_1_8(i_2_1-1) = i_2_1;\n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_8: \" << e.what() << endl; }\n", " double v_1_9; try { v_1_9 = 0;\n", " for(int i_2_0 = 1; i_2_0 <= v_1_8.size(); i_2_0++) {\n", " int i_2_11 = v_1_8[i_2_0-1];\n", " double i_2_12 = v_1_9;\n", " 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)));\n", " double v_2_1; v_2_1 = i_2_12 + v_2_0;\n", " v_1_9 = v_2_1;\n", " \n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_9: \" << e.what() << endl; }\n", " 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; }\n", " double v_1_11; try { v_1_11 = 0;\n", " for(int i_2_0 = 1; i_2_0 <= v_1_8.size(); i_2_0++) {\n", " int i_2_11 = v_1_8[i_2_0-1];\n", " double i_2_12 = v_1_11;\n", " 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)));\n", " double v_2_1; v_2_1 = i_2_12 + v_2_0;\n", " v_1_11 = v_2_1;\n", " \n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_11: \" << e.what() << endl; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " double v_1_15; try { v_1_15 = -v_1_14; } catch(const std::runtime_error& e) { cerr << \"v_1_15: \" << e.what() << endl; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " double v_1_19; try { v_1_19 = min(1, v_1_18); } catch(const std::runtime_error& e) { cerr << \"v_1_19: \" << e.what() << endl; }\n", " 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; }\n", " variant>, tuple>> 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; }\n", " x_cc_0_0 = v_1_21; break;\n", " }\n", " case 1: {\n", " double i_1_0; double i_1_1; Matrix i_1_2; tie(i_1_0, i_1_1, i_1_2) = get<1>(i_0_0);\n", " 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; }\n", " 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; }\n", " Matrix v_1_2; try { v_1_2.resize(v_1_1);\n", " for(int i_2_1 = 1; i_2_1 <= v_1_1; i_2_1++) {\n", " \n", " v_1_2(i_2_1-1) = i_2_1;\n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_2: \" << e.what() << endl; }\n", " double v_1_3; try { v_1_3 = 0;\n", " for(int i_2_0 = 1; i_2_0 <= v_1_2.size(); i_2_0++) {\n", " int i_2_11 = v_1_2[i_2_0-1];\n", " double i_2_12 = v_1_3;\n", " 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)));\n", " double v_2_1; v_2_1 = i_2_12 + v_2_0;\n", " v_1_3 = v_2_1;\n", " \n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_3: \" << e.what() << endl; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " double v_1_8; try { v_1_8 = 0;\n", " for(int i_2_0 = 1; i_2_0 <= v_1_2.size(); i_2_0++) {\n", " int i_2_11 = v_1_2[i_2_0-1];\n", " double i_2_12 = v_1_8;\n", " 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)));\n", " double v_2_1; v_2_1 = i_2_12 + v_2_0;\n", " v_1_8 = v_2_1;\n", " \n", " } } catch(const std::runtime_error& e) { cerr << \"v_1_8: \" << e.what() << endl; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " double v_1_14; try { v_1_14 = -v_1_11; } catch(const std::runtime_error& e) { cerr << \"v_1_14: \" << e.what() << endl; }\n", " 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; }\n", " 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; }\n", " 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; }\n", " double v_1_18; try { v_1_18 = min(1, v_1_17); } catch(const std::runtime_error& e) { cerr << \"v_1_18: \" << e.what() << endl; }\n", " 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; }\n", " variant>, tuple>> 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; }\n", " x_cc_0_0 = v_1_20; break;\n", " }\n", " } } catch(const std::runtime_error& e) { cerr << \"x_cc_0_0: \" << e.what() << endl; }\n", " switch(x_cc_0_0.index()) {\n", " case 0: {\n", " double i_1_1; Matrix i_1_2; tie(i_1_1, i_1_2) = get<0>(x_cc_0_0);\n", " cout << 0 << ' ';\n", " cout << i_1_1 << ' ';\n", " cout << i_1_2.size() << ' ';\n", " for(int i_2_1 = 1; i_2_1 <= i_1_2.size(); i_2_1++) {\n", " cout << i_1_2(i_2_1-1) << ' ';\n", " }\n", " break;\n", " }\n", " case 1: {\n", " double i_1_1; double i_1_2; Matrix i_1_3; tie(i_1_1, i_1_2, i_1_3) = get<1>(x_cc_0_0);\n", " cout << 1 << ' ';\n", " cout << i_1_1 << ' ';\n", " cout << i_1_2 << ' ';\n", " cout << i_1_3.size() << ' ';\n", " for(int i_2_1 = 1; i_2_1 <= i_1_3.size(); i_2_1++) {\n", " cout << i_1_3(i_2_1-1) << ' ';\n", " }\n", " break;\n", " }\n", " } cout << endl;\n", "}\n", "\n", "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" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "compileCC (soccer n `rjmcC` soccerJump)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "-- reduce number of iterations from 5000 to speed this up\n", "samples <- silence' $ iterateLimit 5000 (soccerStep n) (Model1 1 goalData)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Posterior probability in favour of model 1:" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "0.7056588682263547" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "putStrLn \"Posterior probability in favour of model 1:\"\n", "sum [case m of Model1{} -> 1; _ -> 0 | m <- samples] / genericLength samples" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ ":opt svg\n", "import Language.Stochaskell.Plot" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAYAAABQVrO3AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO2dfdAdZZmnr/eFBCIrCCpIgE3YRSsIL0F0S2RhjAHiyAyUqDAlW1uWolWUrKIg67jq5KDDWjq6LlWKhSOM4B8iICKrBcmgBlABvxBP+JTdnJAQEERgrQUJkLN/3H1y+vTpPqc/nu5+us/vqjr1vt399N2//nj66efrvkEIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQggxa+wCfBP4A3BnsDxgNdAF7gcuBOaqFidES1gGrAO2ApuB87H8lCf/KV8K4Zh54DjgGOCXDDPiEuA+YAFYDPwYWFOHQCFawEHAm7BCa3/gAeD1ZM9/ypdCpGQ+Q9odwK3AnyLrVwKPAXcD24GrgJOcqBNi9tgC3Ab0gUexgvBVZM9/ypdCpCRLQZjE/lhzzY5g+WHgQAd2hZh1VgCvxgrGJJLyn/KlECnZtebj96Mr1qxZc/P69etXhdd1u90OsDaS9IKFhYVOJN0G4M0p0qW1p3Szk+7mhYWFVfjDPsDlwLmM1wKrwHXebGQ6gHDaGvVtQO+2sftRJysY7aM4GriFYe3yLOCilLbGMlsRut2uU3tl2PTdXhk2Z1VjSj4KPAjciw2UAdgDuAk4MyZ92vxXJF+C47yZhxrviVcawA8dbdbgomn0LmA/4DCsU/504IYM+/eBjgMdQjSRLwOHAIdio0QXA1cAPwIuTbF/Uv4rmi9BeVPMCFkKwnngGmADcASWac8EngXOBq4E7gHuANY7VSnE7HAccAqWp7YGv1PJnv+UL4VISZY+wh3AuxK23YR9eeZBc5uEGPIjYFHCtqz5r0i+BOVNMSO4aBotisvmlwsc2SnTpu/2yrA5qxqbTt1Noz7cEx80gB86pKEk+tSf2YQQ49Q+MEKIWUGZTQg/0UeqmBna1jQqhBBCNArVCIUYxZXT7SQ7aVHeFKIinDa/BN4InOLapu/2yrA5qxpz4srpdpKdtAR5878+GSxuCm3aFLvHcHsP+vPDtP0e9GMK4YGd/qa4n92TcJq4fcP7T9S0KTndmI7e4P8993zhycl2wc6135twzN74ecTpD+sb3T6uIy5t0jUJ64hLF/0//lrFa4Dh/Y4eJ0rUZtx5w+j1HD03j/KpU1wXhDPnwWQWz7kMmz54zYhhDvgB8DehdWk8y/zPFHam0Q9eUNG/PfsfEgqwYHt/a5CmH/yWDPfZuW8/kmbkt3Tpc/2Y42+K7LtpuM9O6VE9myK2o4VEP3j5xuqIXJa4Qmvr+DWJPa/BsXqjOhLTRbdHr3HSPYruEzmXsXTh/3vj6Xee+2D7fOQc+sNrMKY1cs8G25KuD0B/S9K5lZVPi/oanQM+jE3cncPcRL0XeCKjDbGT6MMXZk7XarbI63T7hBx2Yphbbu4eO8uDFcuH2/qbRpfHOIDRmsGD0N8e2if4m1R7gG3bFg/+jTl+nIb+Jpg7eIqu8LHD+++RnHykZhSn4YAE2wnH3elKb3lMmvD66Pa49NG0U/bZqWt5QppljDB2HsuAhxie82DfAyLpojoi25KuT38TQ+fwyxP+OqdoQbgU+AjwH7DC71KsIPxiBht9bG5Ip6CW1tDppFsnWo0rp9sF7VzAcOrWWkLZdHmKncMv1aWDf7rdbn9hITZNVqIalkN/l5h0U/d/yUt2bHvmmcSxg9HjJGybVGB6w3IH6aOFXh7i7E5aXypFR432Axu7Y001uwHbiooSYoZIcrp9FXAJ8MMp+z+C+RQd5OUDsFphVjsx9CO/TnYTMSwsHN5zYiieh/Ls9Mwz8//GwbGXO7AhaqBoQbgNG6V2d/D/7sB3MtqYQ55lfLJXhs1Z1ZiGspxuZ7UTg8usOcLyMowGLJ2eRDSYUvJp0T6nvYGrgfdjX6GXAbcDX025/1h/2KzEPFtYOJwkkptGR/sIfTiPFqUb214TxwM3Yv1+Az4EfB+r3R2L5bvHMe2XYn2CF2EtMlcDnwTekmDneyl19GOypxA1U844iaJGTwbeh3nHBzgt+P+MlPv3HWjwnOTBL1n6/eIKQiFKJHhuR/oFhaiZct6BRQfLPIKFhNkP+0Jdg/V3ZKH1g2XcDXRJKlRVQAohRF6KFoS/xkak/SJYvhP4eEGbIoakwlSjSUV5qGlUzAYuRo1+Bhvttgx4O9mHaDvrkW+CB5Mm0ASvLU3Q2HxKGywjRC7Kyqdtc7odHSDhq03fcX3OTbgvs3ifp+AyawrhhFLyqQ8FoRDCS1QjFLOBDwWhcpsQQ5bhJvrEgEXAT4Gb8SO/C+EdPmQMtb8IMWQHlh8OwhxqfwA4Cssnl2LTk14IpV8CXIxNWVrAIk6cGNp+JlaovphdirKmmA18KAhnzbNME2iC15YmaMzDFsw5dh94FAuf9CqsgLyV8cFoK4HHMO9O27FJ9ycF25ZiA9guySdFjTXCO7z0LFOUwfjsFs8j7PfLnOKgifat5lDgOqyWNygAVwDfwmqLLzJ0YHFasP1k4D3B8jcwN2vPA/+IeaAZRKmYhjzLCA8p513nQ41QCDFO0egTq4O/t+SXMBf5dfKbEsJjXBSE+2L+Cx/CvMoclXF/5TAxy5QVfeJo4K+BTcB3sVrltWTO82txHX1CCN8oWhDOY80vt2KZ+HVYxsuCcpiYZcqKPnEhViguB96J9Tu+g/RNoyhrilmhaEG4DHgtFm2iD/wZeLKoqLw0wYNJOfT7yb/sNMFrSxM05uQ44BTgbGy051asH3AeuAbYgPn33YyNCH02SHslcA9wB7DejRQ11gi/KCufFu14PBHrhO8BR2KZ8EPA0yn3dxp9wiJfLzjtTE1v002UCZfkHUjj+jrWe1/qsdcCNFhGeEe3u5Ey8mlRp9vz2Nfpx7Cml4uBTwB/n8HGSG4rEo8wSBvNvUXj1UVtxqZbWPDXAXae+Hyhc3YS78+1vZBN13ELxU7mSBGGqYcis4uGU7RkPRwLBHoENkT7JOCD2BDuNJ+TraoR+lYQqkZYn70WkDIe4dxcpDWkR7MLxh7N1t9met3uxuVl5NOifYT3Y30Ur8UKtOODdVnaVNQjL4SX5MmacwfnONBTwd9ezLa4dSWSS3+UngMbLu00ld7oopN7E0vRgvB5rKP+cuD/AK8G/ntRUQVoggeTJtAEry1N0NhwKhssE8QwjXvRlffyKw9Xmp3YCQ9e7E1JO3H7brv1/1JUTDZiz7+1nmXq1uCI9jSNilpZBnwdmw7xIvAV4IvYR+ulwNuAbcAbGPoPXQ1chE2fuAb4FJa39sXmIr4ec792OvCblDpSDpaJNo2ONZWm4RTg+vh9E+09Bbws+L+Hm+bMnr188422HpLrGri208OuyV+x06nCtGbsqcc7H/in8f1GlqPb0uqMYW4O+puG28t7l/ngWUZNo6XiblqFqARXTrcdzPEtpUbYi1m3Lcf+Hx+uKlxz6jmy4wM9+7PzXJ5m/JpH08RsG/sfLPJJzH7h5anXMGJzWvpq7knRUaPCY5JqqL7VXMUIW4IfxDvdXhFJH3a6vYOh0+3fY333pzGc4+sBcwePfuUD5h0nDb1Qje0R0tdEJm0jphY4JX0mHsNq5hn0JLINc6SeRC/mXJ6OrItLEyK8bSzdH3JonmB/Ir3ix0qPDzVCzdoVIp4VWL/7bRPS7I+9oAYeYx4GDsS81TyBeam5P/i7V7bDjzXW9KbsMG17wNhXfsoX7Mh+2xJqIjEa5g5Or23sOGGiNtLYPDKF/TR2AP7thH2eTNAdmdNdqIb1aIa0vQLHoerauQ8FobOm0SZ4MGkCTfDa0gSNBSnqdHswx/crWH/jX7A5vhkYdbq9Zs1bN4c29mD0unW7G5eH5ov2wpb23PPFp4L0/fG5vnMjsRKXLt2+8/+YecEAXHPNg2cMto3aHH+B2vaNy6O2iVzXgZ1ut9uPpBvoPHhwHkuXbqfb3bg8Pt2Qn/3s3qnNvgNtKdK9ED52ZPMplmY493np0u3cddfGJ4N1vXF7w2trNifngQ0b7rsrfIx4jYN7slNf+Li98L7BMWPvb7fb7USflWDdRI15qXsgheYR1sC0QTRNmKPXBI0p+Sg28vp5zEn2Zszp9veBbzPubzQahulo4AvAKqxWeBbmt/SfKTzHt99jfDBF0Kw5GFgRHmAxeKai63fu2x9NE7ctzuaY/U3AIVaARm0mDbiJ3b4SuGtc285LkGRn02iBGzvAZ9D8uwSbYgY7m0R3aoguE2Mn5TXkUJi7b/Q8J5F073bafQrm9g4f46677mblysMm6Ipq3Pk39MyEdceeB6PbRtOXlU9bVSMUooGU5XTbwRzfueXjWTN1s+G09Wm3xzB3cKgWGbd/L7Q+yX6P0WbDWDtJNcMpAsNpnhvaHtmvl9BkG6cjDX/MuV8SH4+umJ9P/ej0RhezDqCpHh8KQiHEEFdOtx3N8X32qelpIPllN61/LnPBmmL/uYOH6yfqejqyPJZm3boH0skY0ovY6MdfgyR9ufvGXAc7eCb424usjy7H4HrkaPm4GjW6CPgJ1lTzFjKFeqm9eVYIn/gRlp/ieFfC+puwGmGU24gdrJGaOWt1/fygGS8tkbR5XnTRkY7ORnEO7AH8X9KdV4bjTyx4s9DLtn6knzVp3ywEBWF0lOfIoKTlkX1SHtfFXE23uKoRnol9ub44LWEMLptGL3Bkp2ybvuP6nJtwX2bxPk8jlDezvMhTp+05tpeSnbWxHSlsX1BPjSXajBq7PrItMc00YmzsrBGGCeWRSr0A9eI1uMNFbWwpcBnweeDTwAmkrxHujEhAY/oJ/Qu3lBV5nBEpCA1kixscEh0MEfcSjA7aiBvYkmQzetwsA0DSLmdh0r7RgSOT9kmy0w+cHQxqS3EDfJJs5TmvwbVNHCT0VzB36+j5JQ0mmqYvfLzo8cPL4RrmmK2EZ8wNRZtG54DPAp/D+iTy2mgUTSnwhChIn1Q1oqL9fKnsFrTlO4Ve8r0SjhdXIyzAtL7Qac2l5dbKixaEq4O/t2CunfIwcvJF4hFWkW5hIeYMGojFUKz++nmcLja+pSjCxJdXz6EtX+jVc9hSCtE0BWEP9323Lu2lpmht7JPY3KTngd0wJ7jrgHeQrnnU6TzCamjOfMEk1DTqNctw43Qb4MPYyNE54F7gvZi3mTRMyJt5m+QmNW8VaeaL2ki7XMR2dNvO6RDh9THnm7WJd1Lzc5HzmXZ8lsPc5tAyCU26jvVNOlZ5FB0scyFwAFaKvxMbpZa2EBzgbLBMEzyYNIEmeG1pgsacuHK6vRT4CHAM8Bpsntl7M2pxPMc3W83Fh3uSXkO5g0e63Y2Xu7KVkv83rmHStaimxl7WM9G2eYTR5jBfbfqO63Nuwn3x5T5vwT4o+8Q73Y66Wws73d7O0Ol2H8vfuwO7YC02WaI8ALX7AfbhnvigAarXEdc0WpWG3oRtpWhwWRD+HGuiyVIbhPozmxC+UsTp9jasxebu4P/dge9kPH4dXp96FR9PjNNj6BZusJzXTg6q7w/2IQxTMDJNhaEQIYo63d4bc7e2EiscL8P8kH41o521wa/sgWw9S7fxchhxtNxPay9uIFs43Vvfup1169LbI2bOWp7jhu1F02YZkBd1UJ3tuFnux7AgsnQbly8sHB7WkNLexlTnS+xzkO5+uKLuARNO5xFW49y5LYNlknDv2FZOtydSltPtm4D3YX2KYHEJT8X6EtMwabDMYA5aqXO7st+TpPmOLjRMGyyT9jj5BsvEX4uyBsvEa+p2N1KvhvLyqQ81QpcnVcYXQ2lfIXWRImCv63Nuwn2p6z5/OfgNKOJ0+36sFvgFbHDMEcG2x4E1wIMZtSW01kzz4+mMjPekFD2+5H8fdLRWQ6tqhNXQ/BphEppW4QXHAzcyGqz2Q1gN8SrgWKzZ83Gs6ehSzJvTRdiAmKuxaU19zNPTmYGNO7EaYtpm1gZObYrisrZSb42w+HGLkugRp9IaYVm0rUYoRNNx6XT7M8EvL03vv+/VLaBF9OoWUCY+TJ+oY2SaEKL1NMIbzTR6dQswWnEtE/GhRiiE8BO11tROuwsgXyhaI1yGuVTbio12O5/smcfZPEK3Xgf6/fhf+2mC15YmaGwBtbbW+HBPatLQi66Y4WtRiYaiBWGSO6gsuMxsTr0OdDrxvxmgCV5bmqCx6dTt7MKHe1KDhtha4IxeizG89CyT5A5KCNF81H8vptGrW4ALXA6WSeMOKo66vzqF8Imk7oZdgG9i0yruDJaZsB7M5WEXm194IerzK0KvbgF+0o4+TFeDZYq4g3IWjzBIG+3Hyx2vbhYZXL/QdXQS78+1vZBN13EL62bQ3XA71rpyM/BjrJC7FAvRdFEofT9h/SAqxWlYQXgjFpVifQYtKjh34uyF33NkRzjExYM+yR3UNJxO2o0LNpuf9k6cT2Ka67W8dt3el3JslqHRAXPA/wK+BvwwWBd1sUbC+jjXayuw0ExpqN3ZhQ/3JJ2G8ieVx+uo3L1ZQ+5HdopexMVYAfgr4HM59q89s00aCTprBWES8jhTG4cC12ExBgctLWkLwoFf0dOC7ScD7yF5Un6UFniWqYq6vKuU6+t1lijaNHoccArwRsxxMJg7qO8VtFspKvCEhxSNPuECZ90WbU5XXhSIaelyRXdoVTpX1P3F58FX5+w1gWZFNcJSKRp9Iml9UlSKc1Lq8qC1pimoZtZ05GJNiHr5MnAIVkhtJnv0iSTCUSkWY1EpbshoQyO6U6FCsOn4UBA6wwfPB+0lv6cdeZbJxKC74WxsCsVWrL9vHrgG2ICFV9qMRZZIWv9sYONK4B7gDrKNGAV5lvFCA/iho80a6m7ucj1qNEfQRjWN5iVtk6kC8zaS2rstfLgnPmjwRUebNfhQI1TTqBB+orwpZgIfok/U/rUlhBBidvGhIKwg+OdsRI2oh6Rrq1GmLUD3UMwEPhSEzjjvvPNuTtqWtx9ww4YOq1bl3LmB9rLYTLqm0fWT7kteXNssQ2MLqDVCvQ/3xAcN4IeONmtw8cW3GvNxuBgbvfYpIhNxJ+C6Qz7BXv4BMZ3OHJ2Ouwql7/Zc2IwZRFPGwIuKnp3KWYb5DT0Mmw/4FeCLWH/+pcDbgG3AG4LtuySsT7JTV97MgzQM8UFHazUUHSwzcOx7BrCAuYI6MaMNdcgLMSQpxufAufapwAuh9EnrfYsVKoS3FC0IVwKPAXcD24GrgJMK2OsUXJ5mjw0bOpmWi+5ftb0ybGbff7K9mARpDEwzmtXmNHt1kRTjcwdwK+Pu1pLWu4wV2pmwnPS/i3STdPisoYx0afZxfVwfNEzT4YSiBeH+WBy0HcHyw8CBGW2EvVdEfctlXSYy2XttdDDHhg0XkGU5Stb9q7bng0brVopOuh+ZgJ/jvk4Nl+Tg2fGOvDE+XdkZ5M1J1y7pfxfpoqSx54OGMtKl2cf1cX3QME2HF5wKXB1aPhnrJ0xLXz/9PPr5xD7AL4C/iaxfAfyS0QC8k9Yn2ZlG3fdCP/2Sfs4pOmr0Ecyf4TxWKzwAqxWmpe6OVyHqJsnp9lXAJQzjEOahiB3lTSFSsgSLfr2AjRr9MZaZhRD5WAx8F/hEwva0NcJpdoQQDjkBGyzzIBac1we3bUI0leOx2uFW4p1uPwo8F6w/c8L6JDtCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEKIgiwD1mGTeDcD5zPusml34KdAD/OUfxGwCPOIEZ4E/LtgXdX6BiwKdN6MTV6uQl9ajUlafNIIsC/wPeAhzPHCUZ5pfC32HA60PBOkq0qjT6wGupjHqAtpn6u1Sc9D0rlnXd80ou84mN1r4ZSDsNiEc1iEigeA10fSzGMvSLBC8SfAKdiL5gFgt5r1DTgLuBJzHzcoCMvWl1ZjkhafNM4D1wPnBuleCuztmcYwi4F7scKxKo2+sAS4j1GXiWtqVeSepOch6dyzrm8i0XfczF2LstyhpYmFtgOLZTjQMXDcXQVpY7UtBd6OOS2uGpfx5MoijcZlWKHy1SDdn4EnPdMY5mgstt/95UvzDtfxRX0k6XlIOves65tG3Dtu5q5FFX5BJ8VCWwz8Bnsg7wNuCNYvATYG686h3Gp2kr454LOY/9TnItuq1DdJ4yQtvmg8BHgCuAIrXK4A9vJMY5i/w3x3vhgsV62xTlzEF20S4ech6dyzrm8SSe+4mbsWZReE+wCXY81i0QjaYF8PR2HNFQcBx2COgo/FHtC3Yg6ET6hB3+rg7y2R9VXqm6YxSYtPGueBI4CvAIcBf8EiIvikccBLsOgp1wbLVWsU1ZHmeWg7Se+4maPMgjBLLLSngX8F/hZrstgcrN8M/AA4sgZ9R2MvxU1YOJs3YS/IuYr0pdGYdK2quoZpND6MDUS5HXgBuA44PKTNB40D3oIN5tkSLFd5HX0gHF8UsscXbQpxz0PSuWdd3ySS3nF/YPauRSkkxUKbx9qOl2ADZZYF6/fEBsuciQ2kGAyieSXwayzyfdX6whzDsCO5Cn1pNSZp8UnjIqz5eyX2EfGl4OeTxgFXYM/ggKo0+sIsxBdNeh6Szj3r+qYSfsfN+rVwRlIstMXY1/ZBWPv8b4FtwbovBdsPB+7BvigeAjq4r7mm0Rcm/JBUoS+txiQtPmkE+9L8LfbleT3wcg81vgyrub4ytG9VGn2i7fFFJ8VpTDr3rOubSPgdB7N9LYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEK7YFbgJ84YeDYkkhMjHvsCvsPyVlvOAc4L/PwnciwVZnnMrTYj68DVUxgvAiVjkaMW1EqJ+FgH/CQvm+l+woMVCtAJfC0KwjLYOeDf6+hTCNQcBG7FgxNcHv72xvPZxLK7cLcDrGAZUXhake3cNeoWYSeaxgvBx4BU1axGiDYSbRg8CnsYCZAOcC1yIFXx3AnsBL8UCKp+D1QjvB3arVrIQ5ZOhr6DfcXvouWn2jsAy4rXAycC/BOv+G/Ab4Atu9QgxczwY/ADWA/8DeBj7AH06WP/9GnQJUSk+N41+EPga8C3gdEzr74CLUf+EaA39TdDvO/xtynDwXSP/D7ogng+t3174FIXwnAw1wqk1OJfsD7wR+AiWKffHmnI2V6hBiAqYO7jGg/97bPDLz7B+vzuA24AvY02hfeAE4Lq6BApRBVmGUVfJ+4ArgWeC5e8D7wS+G2w7ELgduLUWdUK0gweAs4FvAD2sMHwKuBGbvvQ4wyZSIYQQolUcBPwc2KVuIULUjc99hEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGESGYX4JvAH4A7g+UBq4EucD9wITBXtTghhBCibOaB44BjgF8yLAiXAPcBC8Bi4MfAmjoECiGEEFmZz5B2B3Ar8KfI+pXAY8DdwHbgKuAkJ+qEEEKIkslSECaxP9ZcuiNYfhg40IFdIYQQonR2rfn4/eiKNWvW3Lx+/fpV4XXdbrcDrI0kvWBhYaETSbcBeHOKdGntKd3spLt5YWFhFWKA67zZyHQA4bQ16tuA3m1j96NOVjDaR3g0cAvD2uVZwEUpbY1ltiJ0u12n9sqw6bu9MmzOqsaGU/v18OGe+KAB/NDRZg0umkbvAvYDDsMGy5wO3JBh/z7QcaBDCOEW5U0xE2RpGp3HBsIcC+wNbMaqrpcCZwNXArsBVwPr3coUQtSApkGJmSBLQbgDeFfCtpuwGmEelNmEEELUhoum0aK4bH65wJGdMm36bq8Mm7OqsenU3TTqwz3xQQP4oUMaSqJP/ZlNCDFOjkEJ/f7knxB+Uvf0CVDTqBBF2B3rmjgQG8l9LfAx4HnM9eFF2CC2a4BPka2A62Nf4J20O3QSUiatF8IH2tY0KsSssR14B7AceDVwBPA2zPXhxcAZmPvDNwEnZrQ9h/KmmAF8KAiV2YTIzw7MxSFYfp4P1rlwfaiPVDET+FAQOstsgTcCp7i26bu9MmzOqsYKWQz8BngUc4B/A25cH9b6kerDPfFBA/iho80afOgjdMla3Gdc1zZ9t1eGzVnVWBXbgaOAvYBvYxFiXDDSnzjNxdrCwmRjMV5BprrW6na7a9OkS2svazoiz0VdLseAtcG1qPS4vt8PX5gDzgEeAH4PXA+8PMP+3rvJ8t2V1yyecxk2fXAf5YiPAp+nmOtDyDlqtNOJ/+UZNerDPfFBA/iho80aijaNLgU+gn2Bvgb4I/DejDbUDyFEfvYFlgX/7wmcgn2YFnV9WCHx0y0WFg6vW5iYEYo2jfaxwnR3bOj2bsC2oqKEEKnZB3NvuC/wIjYo5ltYc2lR14eZp0/kJW56haZciKooWhBuAy7ERqY9hwXu/U5GGy7nEbbAg8mkJqS5uGvVgnP2wmZTPVbcBxyZsK2I60PQHF/w57nwQUdrNRR90PfGvjTfj41Kuwy4Hfhqyv0V82wsXb+f/HU8WhD6fR6NSxc3QGGWGeTNDDXC+GcX4p/faftN3kcIdxR9yE4G3gecGiyfFvx/Rsr9+w40tAy9FIQX5MibKghFMyk6WOYRzJPFfoGtNcCDGW1osIwQfuI4b8oHqfCTon2EvwYuB34RLN8JfLygzRlBLwDhPYk1uKyG8g98ydxnLkRm6n6QnDaNdrvdjut+Htc2h/aSm5HiSGomatY5u6MJGhtOf9JUwrpHdFbVbOrLc+GDjjZraJWLNcYHSPhos4C92LlWLT/nymyWobHhzLFqVYdOh7HfDOHLc+GDjtZqaJuLtdai8Daiajodtd6L2cCHgrDu5tmSGe3jMH+M6h8U/tPpjGbNVavWsmpVpx4xMeTxXZo1XcwxC9kr4Gu0kvNN4cZZ36UAAA7tSURBVGu0X8dxB+koiboLoRxzlZLpdrv9hYUFp+dU3Ga2vsCslNFX4vo6+nlfyrXXAvo+1wgr7CP04rnwQUebNfjQR+gy1EsTPJg0gSZ4bWmCxkbT6cyxYUOnbhl148tz4YOO1mqo+0vHaY3QT5pXIxQC1QjFDOFDjVAI4SGqEYpZwcVgmX2BS4DXYx7vT8eiZadFX3VCeIjPNUIhXFK0RjgPfAOLOrEMeB2wKaMNuVgrTJLrKo1OFUKIaRQtCJcBr8WiTfSBPwNPFhWVl2DIrfc2XRI32bnopGfX59yE++L7fZ7AMmAdsBXYDJzPsJVlNdAF7sfCpWVqffG/abT8jz9fngsfdLRZQ9GC8BDgCeAKLLNdAeyV0YbLUaNN8GDSBJrgtaUJGqtgB5Z/DgKOBj4AHAUsAS7GIsEsAG8CTsxiuNPpezVvMIzrj78J+PJc+KCjtRqK9hHOY9EnPgbchmW8TwB/n8HGyBdckXiEQVrXk06jNgtNiq2DPJNYQ+fsZFKsa3shm64n7TaNLcEP4FHgAeBVwCLgMSxo9g4scv1JZIhS3+nMeTeJXogyKDpQ5XAsMO8RwPNYRvsgFqcwTROFa6fbHk7cLnf6xCTyDjFvwmT1JmisgUOB67Da35ux2uBpwbaTgfcA70ppy+vpE0m4nlbhy3Phg442ayhaI7wfeBbrJ/wdcHywLksO6tOKeYQamCJqZR8sJNq5wJ9cGPTdxVoa5GLNbbq2ulgrWhA+D5yNZcC9sI759xYVVYBaPZi0yAF2E7y2NEFjVeyBNX1eAvwwWPcIw4DZO4ADgIezGG1ijTBK8ILt5E3X7Xad2subjpRdLWXqmxQCqarrEr0frqi7yu+0abRe6msCTULeN2aCxcC3gV8BnwutXwL8FmsKvR+4EfhC8DcNahoVM4MPnmU0j1CI/BwHnIK1zGwNfqdiXRZnA1cC9wB3kGGgDDRh+oQQbvAhDJMQIj8/wkaIxnETcFhew02sERqT+utVWxTj+FAQNuzBbNqgmCS9eiGIyTRx+sSk7gnfui6EP7SqabQqDyYVTuYtRF6dTfDa0gSNTcfnCfVV4ctz4YOONmvwoSCUZxn/aILXliZobDTqIwT8eS580NFaDT40jbZkHqEQ7aK5fYRCZMOHGqEQwkNUIxSzgqsa4SLgJ8CLwFuwCbxp0aANITxENUIxK7iqEZ6JzV96Mce+LucRNsGDSRNogteWJmgUzceX58IHHa3V4KI2thS4DPg88GngBNLXCHdGJKAxfYT+eZDJirxviBSMVQebNpUiSqcD3e5GqNFXptK598HqgqIvwzksQv0VmN/RfyR7QdiwF7IKQjETNNLF2iT03IskivYRrg7+3oKFfsmDs3iEVaRbWIg5gwYS50DXp+tcQzpv4kf6QhMn1AuRh6JfR5/E4g8+D+wGvAxYB7yDdLVC1QhrQF/GIgWqEYqZoehgmQux8C7LgXdiUerTFoIDGudZpu00wWtLEzQ2HU2f8Oe58EFHmzW0bR5hEzyYNIEmeG1pgsZGIxdrgD/PhQ86WqvBpWeZnzPsM8yCmipqQc64xSyiyBRiHLlYm0GS+jib3vc5w+wCXAq8DdgGvIHhnN7VwEVYAN9rgE8RMzUijrYNllFkCpGEDwWhEKIYfawg/DpW6A1YAlwMnMYwSv2JpAzQ27bBMkIk4UMfocvoE03wYNIEmuC1pQkaq2IHcCvwp8j6lcBjwN3AduAq4KS0RjVYBvDnufBBR2s11N0mLs8yHqHh5Y1nBfAt4GisafRU4AysRghwMvAe4F0pbLVu+sQk9OzPNj40jerhE8JDOp3RrNmm/sI4ut1utOT31fmD0jmm7kJINUKP0Fdx44nWCI8GvgCswppPzwIOBc5JYWsGa4RJKE+0HR9qhEKIcrgL2A84DBssczpWMKaibaNGJ6GR1LNN0cEyyzCXaluBzcD5ZK9lOhss0wQPJk2gCV5bmqCxQuaxqREbgCOwvHgm8CxwNnAlcA9wBylHjIIm1PuED89nmzUULQh3YIXYQVgzzAeAozLacBmPMIfXgX4/+ZfXZuNpgteWJmisih3YAJhXYT5/D8SmUwDchNUIDwE+QTb3hwJI8Y6oAh+ez9ZqKNo0uiX4ATwKPIBlxkYR1/xh6/p9izZR6QMvhBfMUtNoEpqEPxu47CNcAbwac7ydBS87omf3IU8q/DVgYNaYpcEyYrZxVRDuA1wOnMv4pN5pOItHGKTNNAS6LfEFXZB2wECeIdCh++JsSHVg0/UQbRGgGqGYFVx85e8BfB/4NsN+ibQ4jUcYF2w2hYTWTodwRdFpFfnuS7U2y9DYcGZq+kQeqpxu5MPz2WYNRW/iYqwA/BXwuRz7ezCPUAXhNDS/cCbpQ/sn0RdB+aI9FG0aPQ44BXgjNkwb4EPA9wraFULUjGqEYlYoWhD+CFhU0Ia+qIQQQtSGD55lFI9QCA+ZNV+jeShpwJbSyddofjRYphw0WGYm0WCZKWiwTHs0tC0eoYbC+4k8yzQQxSP0Ch+ez9ZqUNOoSEmSdx2NmmsrqhEWYZI3KuUZ3/ChIKzgoZCLtCJMnmivArKtaEJ9GpLfLcmuG4Vv+FAQVlIj1APoHoWuaTeqEU5Gz3l78KGP0BnnnXfeza5tuu4j8d1eGTbLuC+ubZahsenU3UfoQ/+kDxrAj+ezzRpcFISrgS4W+PNCaoxHuH79+je7sBNmw4YLZspeGTbLuC+ubZah0QMK5c264xGW8Wz7oSF7WCcfns82ayjaNLoEuBg4DctsNwInkiH4JxosI0QZuMibwjEK6+QnRWuEK4HHgLuB7cBVwEkF7HUKLjP6dbV27Csr2tQxbTlK1v2rtleGzez7T7YXkyCNgWlGs9qcZq/pFM6bg6bRSfc76X8X6aKkseeDhrzp7JGMvrvCNcW4bdF00W0jxklYTvN/lEnp0trLky7N+swULQj3B/7AMOr1w1h07CyEm0ajc0SyLtPpDH9wQeh/I9rUMW05Stb9q7bng0ar4Cdl4n4fWJtxGabPHyr87LSMwnlz0DQ66X4n/e8iXZQ09nzQkDfd6Pvqgsi7LOm9Fk03ui3EpOc/zf9RJqVLay9PujTrK+dU4OrQ8snANRn27+unn0e/NqG8qV9bfwM6OKJoH+EjwH5YzXIHcAD25ZkWzTUTohyUN0Xb6dQtYMASrCN+AYtN+GPgr2tVJIQA5U0hKuUErEP+QSw4b6vmJgrRYJQ3hRBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEKIFlwDpgK7AZOJ9xl027Az8FesAW4CJgEeYFY2vo97tgXdX6BiwKdN6MTUiuQl9ajUlafNIIsC/wPeAhbHL3UZ5pfC32HA60PBOkq0qjTxSNL+ryOHFpsuTdsjQMiL4bqtYQl6+q1DAHnAM8APweuB54eU4NWTXtAnwTcyx/Z7DsHQcBb8LE749dqNdH0sxjNxKsUPwJcAr2onkA2K1mfQPOAq7EXFQNCsKy9aXVmKTFJ43zWAY5N0j3UmBvzzSGWQzcixWOVWn0hSXAfYy6ZVtT03GS0mS9n2VoGBB9N1SpISlfVanhAGAT8IpAw2XAx3JoyKNpHjgOOAb4JQULwrJcLm0BbsM8hT+KPayviqTZgcVLG+gYOAeugjT6AJYCbwcuqUhXmLQa6ySNxmVYofLVIN2fgSc90xjmaOBP2JforOE6vmiR4ySlcZUvimgAN++GIhpc5asiGvrYe3t3rCDaDdiWQ0MeTTuAW7G8WpgqfA+uAF6NPbxRFgO/wR7o+4AbgvVLgI3BunMo1xN+kr454LOYj8bnItuq1DdJ4yQtvmg8BHgCuAIrXK4A9vJMY5i/w8IVvRgsV62xTlzEF3V1nDRp0tzPMjRMejdUpWFSvqpKwzas2fLu4P/dge/k0JBHk1PKLgj3AS7Hqu9xJfd2rF37oOB3DPA8cCz2gL8VOBNzHly1vtXB31si66vUN01jkhafNM4DRwBfAQ4D/gJ8wjONA16CRWi4NliuWqNIT5r7WRZJ74YqScpXVbI3cDpWgzsA61s/q2INTiizINwDq9JeAvxwStqngX8F/harbm8O1m8GfgAcWYO+o7GX4ibgu1i/xLXY12AV+tJoTLpWVV3DNBofxgai3A68AFwHHB7S5oPGAW/BBh1sCZarvI4+EI5hCNljGLo8zqQ0Wd4tZWhIejdkfZ8W0ZCUr7K2WBTRcCzWJNvDPhqvB/5jxuPn1dQIFmMPSPQLZR5r612CDZRZFqzfExsscyb2lTEYRPNK4NdYdO2q9YU5hmGHeBX60mpM0uKTxkVY8/dKLJN+Kfj5pHHAFdgzOKAqjb5QVQzDpOOE70lSmqT7WaWGMOF3Q5UakvJVlRreAPxvhoXWPwOfyaEhj6YBK3AwWKYsjse+EMJDz0/FTmoL1gy6Avgt1ra8BbuJi7GvmnuwL4CHsCjErmuuafSFCT/sVehLqzFJi08awb6Yf4t9QQ+GWPum8WXYl+0rQ/tWpdEnqophGHec6D2JS5N0P6vUEKZIQVhUQ1y+qlLDHPAPWGvJZqxWuk9ODVk1zWN9+Y9i/bRbGf2IFUIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC+MCuwE2YV/doSCQhhBDCGb6GlHkBOBGLPF1GPDQhhBAC8LcgBIsOvg54N9mjLgshxtkXi1+3Kxbg9CZgr1oVCSEmMo8VhI8Dr6hZixBtYFAQngjcjPKVEIB9Gaak33F76Llp9o4AXgpcC5wM/AvWTHoM8Brg08Dv3WoSomoqz1dLge8AHwX+GKxTvhIzjc9Nox8EvgZ8Czgd03oj8A/AQ8De9UkTorE8B7wTOB+rIYLylRBesj9wF/ASYBHWnLMs2PbB4CeEyEa4j/DD2EfmoP9d+UrMLLvULSCBc4Au8BNgB1Yw/jvgSOD9wBPYiNJH6hIoRAPZA/jPwNeBXwFnA38GVqF8JYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEM3l/wOPCzR4GLLhZAAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "let plot' = layoutToGrid . execEC . plot\n", " hist s vals rng = return . histToPlot $ defaultNormedPlotHist\n", " { _plot_hist_title = s, _plot_hist_values = vals, _plot_hist_range = Just rng }\n", "\n", "toRenderable $ \n", " (plot' (hist \"λ₁\" [real lam | Model1 lam _ <- samples] (2.4,2.7)) `beside`\n", " plot' (line \"lpdf\" [[(i, real $ soccer n `lpdfAuxC` fromConcrete m) | (i,m) <- [1..] `zip` tail samples]]))\n", " `above`\n", " (plot' (hist \"λ₂\" [real lam | Model2 lam kap _ <- samples] (2.4,2.7)) `beside`\n", " plot' (hist \"κ₂\" [real kap | Model2 lam kap _ <- samples] (0,0.08)))" ] } ], "metadata": { "kernelspec": { "display_name": "Haskell", "language": "haskell", "name": "haskell" }, "language_info": { "codemirror_mode": "ihaskell", "file_extension": ".hs", "name": "haskell", "pygments_lexer": "Haskell", "version": "8.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }