{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax, DeriveGeneric, DeriveAnyClass #-} import Language.Stochaskell 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) -- 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 :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)))