{-# LANGUAGE FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax #-}
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
linreg :: Z -> R -> R -> P (RVec,R,R,R,RVec)
linreg n xMu xSigma = do
x <- normals (vector [ xMu | i <- 1...n ]) (vector [ xSigma | i <- 1...n ])
alpha <- normal 0 10
beta <- normal 0 10
sigma <- truncated 0 infinity (cauchy 0 5)
y <- normals (cast alpha + (beta *> x)) (vector [ sigma | i <- 1...n ])
return (x,alpha,beta,sigma,y)
linreg' :: P (Z,RVec,R,R,R,RVec)
linreg' = do
n <- poisson 10
(x,alpha,beta,sigma,y) <- linreg n 0 1
guard $ sigma < 2
return (n,x,alpha,beta,sigma,y)
(nData,xData,alphaTrue,betaTrue,sigmaTrue,yData) <- simulate linreg'
nData
8
:opt svg
import Language.Stochaskell.Plot
toRenderable $ do
plot $ points "data" (list xData `zip` list yData)
plot $ line "truth" [[(x, real alphaTrue + real betaTrue * x) | x <- [-2.5,2.5]]]
posterior = [ (alpha,beta,sigma) | (x,alpha,beta,sigma,y) <- linreg nData 0 1
, x == xData, y == yData ]
samples <- hmcStan 1000 posterior
--- Generating Stan code --- functions { // https://github.com/stan-dev/stan/issues/452 real to_real(real x) { return x; } } data { vector[8] x_stan_0_0; vector[8] x_stan_0_4; } parameters { real x_stan_0_1; real x_stan_0_2; real<lower=0> x_stan_0_3; } model { vector[8] v_0_3; vector[8] v_0_4; vector[8] v_0_6; v_0_3 = x_stan_0_2 * x_stan_0_0; v_0_4 = x_stan_0_1 + v_0_3; for (i_1_1 in 1:8) { v_0_6[i_1_1] = x_stan_0_3; } x_stan_0_1 ~ normal(0, 10); x_stan_0_2 ~ normal(0, 10); x_stan_0_3 ~ cauchy(0, 5) T[0,]; x_stan_0_4 ~ normal(v_0_4, v_0_6); } make -C /home/jovyan/stochaskell/cmdstan /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca make[1]: Entering directory '/home/jovyan/stochaskell/cmdstan' --- Translating Stan model to C++ code --- bin/stanc --o=/home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.hpp /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.stan Model name=model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca_model Input file=/home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.stan Output file=/home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.hpp --- Compiling, linking C++ code --- g++ -std=c++1y -pthread -Wno-sign-compare -O3 -I src -I stan/src -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.69.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -c -x c++ -o /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.o /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.hpp g++ -std=c++1y -pthread -Wno-sign-compare -O3 -I src -I stan/src -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.69.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION src/cmdstan/main.o stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_idas.a /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca.o -o /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca make[1]: Leaving directory '/home/jovyan/stochaskell/cmdstan' --- Sampling Stan model --- /home/jovyan/stochaskell/cache/stan/model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca method=sample num_samples=1000 num_warmup=1000 save_warmup=0 thin=1 adapt engaged=1 algorithm=hmc engine=nuts max_depth=10 metric=diag_e stepsize=1.0 stepsize_jitter=0.0 data file=/tmp/stan-8bebf31206ad0499/stan.data init=0 output file=/tmp/stan-8bebf31206ad0499/stan.csv method = sample (Default) sample num_samples = 1000 (Default) num_warmup = 1000 (Default) save_warmup = 0 (Default) thin = 1 (Default) adapt engaged = 1 (Default) gamma = 0.050000000000000003 (Default) delta = 0.80000000000000004 (Default) kappa = 0.75 (Default) t0 = 10 (Default) init_buffer = 75 (Default) term_buffer = 50 (Default) window = 25 (Default) algorithm = hmc (Default) hmc engine = nuts (Default) nuts max_depth = 10 (Default) metric = diag_e (Default) metric_file = (Default) stepsize = 1 (Default) stepsize_jitter = 0 (Default) id = 0 (Default) data file = /tmp/stan-8bebf31206ad0499/stan.data init = 0 random seed = 1482352849 output file = /tmp/stan-8bebf31206ad0499/stan.csv diagnostic_file = (Default) refresh = 100 (Default) Gradient evaluation took 1.4e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds. Adjust your expectations accordingly! Iteration: 1 / 2000 [ 0%] (Warmup) Iteration: 100 / 2000 [ 5%] (Warmup) Iteration: 200 / 2000 [ 10%] (Warmup) Iteration: 300 / 2000 [ 15%] (Warmup) Iteration: 400 / 2000 [ 20%] (Warmup) Iteration: 500 / 2000 [ 25%] (Warmup) Iteration: 600 / 2000 [ 30%] (Warmup) Iteration: 700 / 2000 [ 35%] (Warmup) Iteration: 800 / 2000 [ 40%] (Warmup) Iteration: 900 / 2000 [ 45%] (Warmup) Iteration: 1000 / 2000 [ 50%] (Warmup) Iteration: 1001 / 2000 [ 50%] (Sampling) Iteration: 1100 / 2000 [ 55%] (Sampling) Iteration: 1200 / 2000 [ 60%] (Sampling) Iteration: 1300 / 2000 [ 65%] (Sampling) Iteration: 1400 / 2000 [ 70%] (Sampling) Iteration: 1500 / 2000 [ 75%] (Sampling) Iteration: 1600 / 2000 [ 80%] (Sampling) Iteration: 1700 / 2000 [ 85%] (Sampling) Iteration: 1800 / 2000 [ 90%] (Sampling) Iteration: 1900 / 2000 [ 95%] (Sampling) Iteration: 2000 / 2000 [100%] (Sampling) Elapsed Time: 0.023748 seconds (Warm-up) 0.030511 seconds (Sampling) 0.054259 seconds (Total) Stan took 0.056999469s # stan_version_major = 2 # stan_version_minor = 20 # stan_version_patch = 0 # model = model_10a646b17ebb76a1b1f1fa35d92ec91da963e6ca_model # method = sample (Default) # sample # num_samples = 1000 (Default) # num_warmup = 1000 (Default) # save_warmup = 0 (Default) # thin = 1 (Default) # adapt # engaged = 1 (Default) # gamma = 0.050000000000000003 (Default) # delta = 0.80000000000000004 (Default) # kappa = 0.75 (Default) # t0 = 10 (Default) # init_buffer = 75 (Default) # term_buffer = 50 (Default) # window = 25 (Default) # algorithm = hmc (Default) # hmc # engine = nuts (Default) # nuts # max_depth = 10 (Default) # metric = diag_e (Default) # metric_file = (Default) # stepsize = 1 (Default) # stepsize_jitter = 0 (Default) # id = 0 (Default) # data # file = /tmp/stan-8bebf31206ad0499/stan.data # init = 0 # random # seed = 1482352849 # output # file = /tmp/stan-8bebf31206ad0499/stan.csv # diagnostic_file = (Default) # refresh = 100 (Default) # Adaptation terminated # Step size = 0.515356 # Diagonal elements of inverse mass matrix: # 0.0259429, 0.0251066, 0.119303 # # Elapsed Time: 0.023748 seconds (Warm-up) # 0.030511 seconds (Sampling) # 0.054259 seconds (Total) # Extracting: x_stan_0_1, x_stan_0_2, x_stan_0_3 --- Removing temporary files ---
toRenderable $ do
plot $ points "data" (list xData `zip` list yData)
plot $ line "truth" [[(x, real alphaTrue + real betaTrue * x) | x <- [-2.5,2.5]]]
setColors [black `withOpacity` 0.1]
plot $ line "posterior" [[(x, real a + real b * x) | x <- [-2.5,2.5]] | (i,(a,b,s)) <- [0..] `zip` samples, i `mod` 100 == 0]
measurementError :: Z -> R -> P (R,R,RVec,RVec,R,R,R,RVec)
measurementError n tau = do
xMu <- uniform (-10) 10
xSigma <- uniform 0 10
(x,alpha,beta,sigma,y) <- linreg n xMu xSigma
xMeas <- normals x (vector [ tau | i <- 1...n ])
return (xMu,xSigma,x,xMeas,alpha,beta,sigma,y)
(nTrue,xMu,xSigma,xTrue,xMeasData,alphaTrue,betaTrue,yData) <- simulate $ do
n <- poisson 10
(xMu,xSigma,x,xMeas,alpha,beta,sigma,y) <- measurementError n 1
guard $ sigma < 2
return (n,xMu,xSigma,x,xMeas,alpha,beta,y)
rejecting OOB sample -12.046149701419338 rejecting OOB sample -6.697030555008025
toRenderable $ do
plot $ points "data" (list xMeasData `zip` list yData)
plot $ points "truth" (list xTrue `zip` list yData)
plot $ line "truth" [[(x, real alphaTrue + real betaTrue * x) | x <- [real xMu - 2.5 * real xSigma, real xMu + 2.5 * real xSigma]]]
samples <- hmcStanInit 1000 [ (xMu,xSigma,x,alpha,beta,sigma)
| (xMu,xSigma,x,xMeas,alpha,beta,sigma,y) <- measurementError nTrue 1
, xMeas == xMeasData, y == yData ]
(0,1,xMeasData,0,0,1)
--- Generating Stan code --- functions { // https://github.com/stan-dev/stan/issues/452 real to_real(real x) { return x; } } data { vector[12] x_stan_0_6; vector[12] x_stan_0_7; } parameters { real x_stan_0_0; real x_stan_0_1; vector[12] x_stan_0_2; real x_stan_0_3; real x_stan_0_4; real<lower=0> x_stan_0_5; } model { vector[12] v_0_1; vector[12] v_0_4; vector[12] v_0_5; vector[12] v_0_7; vector[12] v_0_9; vector[12] v_0_11; for (i_1_1 in 1:12) { v_0_1[i_1_1] = 1; } v_0_4 = x_stan_0_4 * x_stan_0_2; v_0_5 = x_stan_0_3 + v_0_4; for (i_1_1 in 1:12) { v_0_7[i_1_1] = x_stan_0_5; } for (i_1_1 in 1:12) { v_0_9[i_1_1] = x_stan_0_0; } for (i_1_1 in 1:12) { v_0_11[i_1_1] = x_stan_0_1; } x_stan_0_0 ~ uniform(-10, 10); x_stan_0_1 ~ uniform(0, 10); x_stan_0_2 ~ normal(v_0_9, v_0_11); x_stan_0_3 ~ normal(0, 10); x_stan_0_4 ~ normal(0, 10); x_stan_0_5 ~ cauchy(0, 5) T[0,]; x_stan_0_6 ~ normal(v_0_5, v_0_7); x_stan_0_7 ~ normal(x_stan_0_2, v_0_1); } make -C /home/jovyan/stochaskell/cmdstan /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595 make[1]: Entering directory '/home/jovyan/stochaskell/cmdstan' --- Translating Stan model to C++ code --- bin/stanc --o=/home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.hpp /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.stan Model name=model_a7a3b527bcc956e0c332c9fd2faa65b486c94595_model Input file=/home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.stan Output file=/home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.hpp --- Compiling, linking C++ code --- g++ -std=c++1y -pthread -Wno-sign-compare -O3 -I src -I stan/src -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.69.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -c -x c++ -o /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.o /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.hpp g++ -std=c++1y -pthread -Wno-sign-compare -O3 -I src -I stan/src -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.3.3 -I stan/lib/stan_math/lib/boost_1.69.0 -I stan/lib/stan_math/lib/sundials_4.1.0/include -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION src/cmdstan/main.o stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_4.1.0/lib/libsundials_idas.a /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595.o -o /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595 make[1]: Leaving directory '/home/jovyan/stochaskell/cmdstan' --- Sampling Stan model --- /home/jovyan/stochaskell/cache/stan/model_a7a3b527bcc956e0c332c9fd2faa65b486c94595 method=sample num_samples=1000 num_warmup=1000 save_warmup=0 thin=1 adapt engaged=1 algorithm=hmc engine=nuts max_depth=10 metric=diag_e stepsize=1.0 stepsize_jitter=0.0 data file=/tmp/stan-350dc34169e2e83f/stan.data init=/tmp/stan-350dc34169e2e83f/stan.init output file=/tmp/stan-350dc34169e2e83f/stan.csv method = sample (Default) sample num_samples = 1000 (Default) num_warmup = 1000 (Default) save_warmup = 0 (Default) thin = 1 (Default) adapt engaged = 1 (Default) gamma = 0.050000000000000003 (Default) delta = 0.80000000000000004 (Default) kappa = 0.75 (Default) t0 = 10 (Default) init_buffer = 75 (Default) term_buffer = 50 (Default) window = 25 (Default) algorithm = hmc (Default) hmc engine = nuts (Default) nuts max_depth = 10 (Default) metric = diag_e (Default) metric_file = (Default) stepsize = 1 (Default) stepsize_jitter = 0 (Default) id = 0 (Default) data file = /tmp/stan-350dc34169e2e83f/stan.data init = /tmp/stan-350dc34169e2e83f/stan.init random seed = 1482372444 output file = /tmp/stan-350dc34169e2e83f/stan.csv diagnostic_file = (Default) refresh = 100 (Default) Gradient evaluation took 4.3e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds. Adjust your expectations accordingly! Iteration: 1 / 2000 [ 0%] (Warmup) Iteration: 100 / 2000 [ 5%] (Warmup) Iteration: 200 / 2000 [ 10%] (Warmup) Iteration: 300 / 2000 [ 15%] (Warmup) Iteration: 400 / 2000 [ 20%] (Warmup) Iteration: 500 / 2000 [ 25%] (Warmup) Iteration: 600 / 2000 [ 30%] (Warmup) Iteration: 700 / 2000 [ 35%] (Warmup) Iteration: 800 / 2000 [ 40%] (Warmup) Iteration: 900 / 2000 [ 45%] (Warmup) Iteration: 1000 / 2000 [ 50%] (Warmup) Iteration: 1001 / 2000 [ 50%] (Sampling) Iteration: 1100 / 2000 [ 55%] (Sampling) Iteration: 1200 / 2000 [ 60%] (Sampling) Iteration: 1300 / 2000 [ 65%] (Sampling) Iteration: 1400 / 2000 [ 70%] (Sampling) Iteration: 1500 / 2000 [ 75%] (Sampling) Iteration: 1600 / 2000 [ 80%] (Sampling) Iteration: 1700 / 2000 [ 85%] (Sampling) Iteration: 1800 / 2000 [ 90%] (Sampling) Iteration: 1900 / 2000 [ 95%] (Sampling) Iteration: 2000 / 2000 [100%] (Sampling) Elapsed Time: 0.246155 seconds (Warm-up) 0.466038 seconds (Sampling) 0.712193 seconds (Total) Stan took 0.720883952s # stan_version_major = 2 # stan_version_minor = 20 # stan_version_patch = 0 # model = model_a7a3b527bcc956e0c332c9fd2faa65b486c94595_model # method = sample (Default) # sample # num_samples = 1000 (Default) # num_warmup = 1000 (Default) # save_warmup = 0 (Default) # thin = 1 (Default) # adapt # engaged = 1 (Default) # gamma = 0.050000000000000003 (Default) # delta = 0.80000000000000004 (Default) # kappa = 0.75 (Default) # t0 = 10 (Default) # init_buffer = 75 (Default) # term_buffer = 50 (Default) # window = 25 (Default) # algorithm = hmc (Default) # hmc # engine = nuts (Default) # nuts # max_depth = 10 (Default) # metric = diag_e (Default) # metric_file = (Default) # stepsize = 1 (Default) # stepsize_jitter = 0 (Default) # id = 0 (Default) # data # file = /tmp/stan-350dc34169e2e83f/stan.data # init = /tmp/stan-350dc34169e2e83f/stan.init # random # seed = 1482372444 # output # file = /tmp/stan-350dc34169e2e83f/stan.csv # diagnostic_file = (Default) # refresh = 100 (Default) # Adaptation terminated # Step size = 0.0422756 # Diagonal elements of inverse mass matrix: # 2.15764, 1.55568, 0.173876, 0.218672, 0.312825, 0.279269, 0.201198, 0.232343, 0.181752, 0.279712, 0.228128, 0.239036, 0.246002, 0.20854, 56.8505, 0.8343, 0.313135 # # Elapsed Time: 0.246155 seconds (Warm-up) # 0.466038 seconds (Sampling) # 0.712193 seconds (Total) # Extracting: x_stan_0_0, x_stan_0_1, x_stan_0_2, x_stan_0_3, x_stan_0_4, x_stan_0_5 --- Removing temporary files ---
let endpoints = [real xMu - 2.5 * real xSigma, real xMu + 2.5 * real xSigma]
toRenderable $ do
plot $ points "data" (list xMeasData `zip` list yData)
plot $ points "truth" (list xTrue `zip` list yData)
plot $ line "truth" [[(x, real alphaTrue + real betaTrue * x) | x <- endpoints]]
setColors [black `withOpacity` 0.1]
plot $ line "posterior" [[(x, real a + real b * x) | x <- endpoints] | (i,(_,_,_,a,b,_)) <- [0..] `zip` samples, i `mod` 100 == 0]