using StatisticalRethinking
gr(size=(500,800))

Plots.GRBackend()

ProjDir = rel_path("..", "chapters", "03")
cd(ProjDir)

binomialstanmodel = "
// Inferring a Rate
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
"

"\n// Inferring a Rate\ndata {\n  int N;\n  int<lower=0> k[N];\n  int<lower=1> n[N];\n}\nparameters {\n  real<lower=0,upper=1> theta;\n  real<lower=0,upper=1> thetaprior;\n}\nmodel {\n  // Prior Distribution for Rate Theta\n  theta ~ beta(1, 1);\n  thetaprior ~ beta(1, 1);\n\n  // Observed Counts\n  k ~ binomial(n, theta);\n}\n"

stanmodel = Stanmodel(name="binomial", monitors = ["theta"], model=binomialstanmodel,
  output_format=:mcmcchain)

N2 = 4^2
d = Binomial(9, 0.66)
n2 = Int.(9 * ones(Int, N2))
k2 = rand(d, N2)

=====> /home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03


File /home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.stan will be updated.

16-element Array{Int64,1}:
 6
 6
 5
 6
 7
 7
 9
 6
 5
 6
 5
 7
 4
 5
 5
 8

binomialdata = [
  Dict("N" => length(n2), "n" => n2, "k" => k2)
]

1-element Array{Dict{String,Any},1}:
 Dict("N"=>16,"k"=>[6, 6, 5, 6, 7, 7, 9, 6, 5, 6, 5, 7, 4, 5, 5, 8],"n"=>[9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9])

rc, chn, cnames = stan(stanmodel, binomialdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME)

describe(chn)

--- Translating Stan model to C++ code ---
bin/stanc  /home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.stan --o=/home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.hpp
Model name=binomial_model
Input file=/home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.stan
Output file=/home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.hpp

--- Linking C++ model ---
g++ -Wall -I . -isystem stan/lib/stan_math/lib/eigen_3.3.3 -isystem stan/lib/stan_math/lib/boost_1.66.0 -isystem stan/lib/stan_math/lib/sundials_3.1.0/include -std=c++1y -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -Wno-unused-function -Wno-uninitialized -I src -isystem stan/src -isystem stan/lib/stan_math/ -DFUSION_MAX_VECTOR_SIZE=12 -Wno-unused-local-typedefs -DEIGEN_NO_DEBUG -DNO_FPRINTF_OUTPUT -pipe   src/cmdstan/main.cpp  -O3 -o /home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial -include /home/travis/build/StanJulia/StatisticalRethinking.jl/docs/build/03/tmp/binomial.hpp stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_cvodes.a stan/lib/stan_math/lib/sundials_3.1.0/lib/libsundials_idas.a

Length of data array is not equal to nchains,
all chains will use the first data dictionary.

Calling /home/travis/cmdstan/bin/stansummary to infer across chains.

Inference for Stan model: binomial_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.

Warmup took (0.026, 0.026, 0.025, 0.026) seconds, 0.10 seconds total
Sampling took (0.037, 0.040, 0.039, 0.030) seconds, 0.15 seconds total

                Mean     MCSE  StdDev     5%   50%   95%    N_Eff  N_Eff/s    R_hat
lp__             -95  3.1e-02     1.1    -97   -95   -94  1.4e+03  9.5e+03  1.0e+00
accept_stat__   0.91  1.8e-03    0.12   0.66  0.96   1.0  4.4e+03  3.0e+04  1.0e+00
stepsize__      0.85  4.7e-02   0.066   0.77  0.87  0.95  2.0e+00  1.4e+01  3.6e+13
treedepth__      2.0  7.1e-02    0.71    1.0   2.0   3.0  1.0e+02  6.8e+02  1.0e+00
n_leapfrog__     4.8  2.8e-01     5.9    1.0   3.0    15  4.3e+02  3.0e+03  1.0e+00
divergent__     0.00     -nan    0.00   0.00  0.00  0.00     -nan     -nan     -nan
energy__          96  4.2e-02     1.5     94    96    99  1.3e+03  8.8e+03  1.0e+00
theta           0.67  7.6e-04   0.039   0.61  0.67  0.73  2.7e+03  1.8e+04  1.0e+00
thetaprior      0.49  5.2e-03    0.29  0.041  0.49  0.95  3.2e+03  2.2e+04  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).

Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
         Mean        SD       Naive SE       MCSE      ESS
theta 0.66927933 0.03934901 0.0006221625 0.0007931537 1000

Quantiles:
         2.5%      25.0%    50.0%    75.0%     97.5%
theta 0.59246722 0.641689 0.669947 0.6970875 0.744323

if rc == 0
  plot(chn)
end