Run Stan models (normal mixture)

.travis.yml

@@ -37,6 +37,7 @@ script:
              if ENV["GROUP"] != "Bench"
                 Pkg.test("Turing"; coverage=true)
              else
+                include(Pkg.dir("Turing")*"/benchmarks/installstan.jl");
                 include(Pkg.dir("Turing")*"/benchmarks/benchmark.jl")
              end'
 after_success:

benchmarks/benchmark.jl

@@ -6,9 +6,12 @@ CONFIG = Dict(
     #"naive-bayes",
     #"normal-loc",
     "simple-normal-mixture",
+    "simple-normal-mixture-stan",
     "simplegauss",
+    "simplegauss-stan",
     "gauss",
     "bernoulli",
+    "bernoulli-stan",
     "gdemo-geweke"
     #"negative-binomial"
   ],

benchmarks/bernoulli-stan.data.jl

@@ -0,0 +1,6 @@
+const berstandata = [
+  Dict(
+  "N" => 10,
+  "y" => [0,1,0,0,0,0,0,0,0,1]
+  )
+]

benchmarks/bernoulli-stan.model.jl

@@ -0,0 +1,14 @@
+const berstanmodel = "
+data {
+  int<lower=0> N;
+  int<lower=0,upper=1> y[N];
+}
+parameters {
+  real<lower=0,upper=1> theta;
+}
+model {
+  theta ~ beta(1,1);
+  for (n in 1:N)
+    y[n] ~ bernoulli(theta);
+}
+"

benchmarks/bernoulli-stan.run.jl

@@ -0,0 +1,15 @@
+using Mamba, Stan
+
+include("bernoulli-stan.data.jl")
+include("bernoulli-stan.model.jl")
+
+berstan = Stanmodel(name="bernoulli", model=berstanmodel, nchains=1);
+
+ber_stan_sim = stan(berstan, berstandata, CmdStanDir=CMDSTAN_HOME)
+describe(ber_stan_sim)
+
+theta_stan = ber_stan_sim[1:1000, ["theta"], :].value[:]
+theta_turing = ber_sim[:theta]
+
+println("The difference between the means from Stan and Turing is:")
+println("  E(theta_stan) - E(theta_turing) = $(mean(theta_stan) - mean(theta_turing))")

benchmarks/bernoulli.run.jl

@@ -5,6 +5,6 @@ using Base.Test
 include("bernoulli.data.jl")
 include("bernoulli.model.jl")
 
-sim1 = sample(bermodel(berdata), HMC(2000, 0.25, 5))
+ber_sim = sample(bermodel(berdata), Turing.HMC(1000, 0.25, 5))
 
-describe(sim1)
+describe(ber_sim)

benchmarks/installstan.jl

@@ -0,0 +1,5 @@
+run(`wget https://github.com/stan-dev/cmdstan/releases/download/v2.14.0/cmdstan-2.14.0.zip`)
+run(`unzip -qq cmdstan-2.14.0.zip`)
+CMDSTAN_HOME = pwd() * "/cmdstan-2.14.0"
+println("CMDStan is installed at path: $CMDSTAN_HOME")
+Pkg.add("Stan")

benchmarks/simple-normal-mixture-stan.model.jl

@@ -0,0 +1,29 @@
+# normal mixture, unknown proportion and means, known variance
+# p(y|mu,theta) = theta * Normal(y|mu[1],1) + (1-theta) * Normal(y|mu[2],1);
+
+const simplenormalmixturemodel = "
+data {
+  int<lower=0>  N;
+  real y[N];
+}
+parameters {
+  real<lower=0,upper=1> theta;
+  real mu[2];
+}
+transformed parameters {
+  real log_theta;
+  real log_one_minus_theta;
+
+  log_theta <- log(theta);
+  log_one_minus_theta <- log(1.0 - theta);
+}
+model {
+  theta ~ uniform(0,1); // equivalently, ~ beta(1,1);
+  for (k in 1:2)
+    mu[k] ~ normal(0,10);
+  for (n in 1:N)
+    increment_log_prob(
+      log_sum_exp(log_theta + normal_log(y[n],mu[1],1.0),
+                  log_one_minus_theta + normal_log(y[n],mu[2],1.0)));
+}
+"

benchmarks/simple-normal-mixture-stan.run.jl

@@ -0,0 +1,15 @@
+using Mamba, Stan
+
+include("simple-normal-mixture-stan.model.jl")
+include("simple-normal-mixture-stan.data.jl")
+
+simplenormalmixturestan = Stanmodel(name="normalmixture", model=simplenormalmixturemodel, nchains=1);
+
+nm_stan_sim = stan(simplenormalmixturestan, simplenormalmixturestandata, CmdStanDir=CMDSTAN_HOME)
+# describe(nm_stan_sim)
+
+mu_1 = nm_stan_sim[1:1000, ["mu.1"], :].value[:]
+mu_2 = nm_stan_sim[1:1000, ["mu.2"], :].value[:]
+println("Result from Stan for Simple Normal Mixture model")
+println("means from 1000 iterations:")
+println("[$(mean(mu_1)),$(mean(mu_2))]")

benchmarks/simple-normal-mixture.run.jl

@@ -8,5 +8,6 @@ include("simple-normal-mixture.model.jl")
 # NOTE: I only run a sub-set of the data as running the whole is quite slow
 nmchain = sample(nmmodel(y[1:100]), Gibbs(250, PG(20, 1, :k), HMC(1, 0.2, 3, :mu, :theta)))
 # describe(nmchain)
-println("means:")
+println("Result from Turing for Simple Normal Mixture model")
+println("means, 250 iteration from Gibbs with PG(20, 1, :k), HMC(1, 0.2, 3, :mu, :theta):")
 println(mean([[Float64(n) for n in ns] for ns in nmchain[:mu]]))

benchmarks/simplegauss-stan.data.jl

@@ -0,0 +1,6 @@
+const simplegaussstandata = [
+  Dict(
+  "N" => 2,
+  "x" => [1.5, 2.0]
+  )
+]

benchmarks/simplegauss-stan.model.jl

@@ -0,0 +1,16 @@
+
+const simplegaussstanmodel = "
+data {
+  int<lower=0>  N;
+  real x[N];
+}
+parameters {
+  real<lower=0> s;
+  real m;
+}
+model {
+  s ~ inv_gamma(2,3);
+  m ~ normal(0,sqrt(s));
+    x ~ normal(m,sqrt(s));
+}
+"

benchmarks/simplegauss-stan.run.jl

@@ -0,0 +1,42 @@
+using Mamba, Stan
+
+include("simplegauss-stan.data.jl")
+include("simplegauss-stan.model.jl")
+
+simplegaussstan = Stanmodel(name="simplegauss", model=simplegaussstanmodel, nchains=1);
+
+simple_gauss_stan_sim = stan(simplegaussstan, simplegaussstandata, CmdStanDir=CMDSTAN_HOME)
+describe(simple_gauss_stan_sim)
+
+s_stan = simple_gauss_stan_sim[1:1000, ["s"], :].value[:]
+m_stan = simple_gauss_stan_sim[1:1000, ["m"], :].value[:]
+
+println("Correctness check for Stan:")
+
+print("  1. s ≈ 49/24 (ϵ = 0.15)")
+ans1 = abs(mean(s_stan) - 49/24) <= 0.15
+if ans1
+  print_with_color(:green, " ✓\n")
+else
+  print_with_color(:red, " X\n")
+  print_with_color(:red, "    s = $(mean(simple_gauss_sim[:s])), diff = $(abs(mean(simple_gauss_sim[:s]) - 49/24))\n")
+end
+
+print("  2. m ≈ 7/6 (ϵ = 0.15)")
+ans2 = abs(mean(m_stan) - 7/6) <= 0.15
+if ans2
+  print_with_color(:green, "   ✓\n")
+else
+  print_with_color(:red, "   X\n")
+  print_with_color(:red, "     m = $(mean(simple_gauss_sim[:m])), diff = $(abs(mean(simple_gauss_sim[:m]) - 7/6))\n")
+end
+
+
+
+s_turing = simple_gauss_sim[:s]
+m_turing = simple_gauss_sim[:m]
+
+
+println("The difference between the means from Stan and Turing is:")
+println("  E(s_stan) - E(s_turing) = $(mean(s_stan) - mean(s_turing))")
+println("  E(m_stan) - E(m_turing) = $(mean(m_stan) - mean(m_turing))")

benchmarks/simplegauss.run.jl

@@ -4,24 +4,26 @@ using Base.Test
 
 include("simplegauss.data.jl")
 include("simplegauss.model.jl")
-alg = Gibbs(2000, PG(30, 3, :s), HMC(2, 0.1, 3, :m))
-sim = sample(simplegaussmodel(simplegaussdata), alg)
-describe(sim)
+alg = Gibbs(1000, PG(100, 3, :s), HMC(2, 0.1, 3, :m))
+simple_gauss_sim = sample(simplegaussmodel(simplegaussdata), alg)
+describe(simple_gauss_sim)
+
+println("Correctness check for Turing:")
 
 print("  1. s ≈ 49/24 (ϵ = 0.15)")
-ans1 = abs(mean(sim[:s]) - 49/24) <= 0.15
+ans1 = abs(mean(simple_gauss_sim[:s]) - 49/24) <= 0.15
 if ans1
   print_with_color(:green, " ✓\n")
 else
   print_with_color(:red, " X\n")
-  print_with_color(:red, "    s = $(mean(sim[:s])), diff = $(abs(mean(sim[:s]) - 49/24))\n")
+  print_with_color(:red, "    s = $(mean(simple_gauss_sim[:s])), diff = $(abs(mean(simple_gauss_sim[:s]) - 49/24))\n")
 end
 
 print("  2. m ≈ 7/6 (ϵ = 0.15)")
-ans2 = abs(mean(sim[:m]) - 7/6) <= 0.15
+ans2 = abs(mean(simple_gauss_sim[:m]) - 7/6) <= 0.15
 if ans2
   print_with_color(:green, "   ✓\n")
 else
   print_with_color(:red, "   X\n")
-  print_with_color(:red, "     m = $(mean(sim[:m])), diff = $(abs(mean(sim[:m]) - 7/6))\n")
+  print_with_color(:red, "     m = $(mean(simple_gauss_sim[:m])), diff = $(abs(mean(simple_gauss_sim[:m]) - 7/6))\n")
 end