search

stan-dev / example-models

Example models for Stan
http://mc-stan.org/
775 stars 477 forks source link

https://github.com/stan-dev/example-models/blob/master/Bayesian_Cognitive_Modeling/ParameterEstimation/Binomial/Rate_5_Stan.R #13

Closed drewyallop closed 9 years ago

drewyallop commented 9 years ago

Warnings and errors. See below

rstan (Version 2.5.0, packaged: 2014-10-22 14:19:22 UTC, GitRev: e52c66f42e81)

to be passed on to Stan

data <- read_rdump("Rate_5.data.R")

myinits <- list(

  • list(theta=.5))

    parameters to be monitored:

parameters <- c("theta", "postpredk1", "postpredk2")

The following command calls Stan with specific options.

For a detailed description type "?stan".

samples <- stan(file="Rate_5_model.stan",

  • data=data,
  • init=myinits, # If not specified, gives random inits
  • pars=parameters,
  • iter=10000,
  • chains=1,
  • thin=1,

  • warmup = 100, # Stands for burn-in; Default = iter/2

  • seed = 123 # Setting seed; Default is random seed

  • )

TRANSLATING MODEL 'Rate_5_model' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'Rate_5_model' NOW. In file included from file1cdfd6e9e25.cpp:407: In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include/rstan/rstaninc.hpp:3: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include/rstan/stan_fit.hpp:732:15: warning: unused variable 'return_code' [-Wunused-variable] int return_code = stan::common::do_bfgs_optimize(model, lbfgs, base_rng, ^ /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include/rstan/stan_fit.hpp:796:15: warning: unused variable 'return_code' [-Wunused-variable] int return_code = stan::common::do_bfgs_optimize(model, bfgs, base_rng, ^ /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include/rstan/stan_fit.hpp:616:14: warning: unused variable 'init_log_prob' [-Wunused-variable] double init_log_prob; ^ In file included from file1cdfd6e9e25.cpp:8: In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:16: In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/agrad/rev.hpp:5: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/agrad/rev/chainable.hpp:87:17: warning: 'static' function 'set_zero_all_adjoints' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration] static void set_zero_all_adjoints() { ^ In file included from file1cdfd6e9e25.cpp:8: In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:23: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/io/dump.hpp:26:14: warning: function 'product' is not needed and will not be emitted [-Wunneeded-internal-declaration] size_t product(std::vector dims) { ^ 5 warnings generated.

SAMPLING FOR MODEL 'Rate_5_model' NOW (CHAIN 1).

Iteration: 1 / 10000 0% Iteration: 1000 / 10000 10% Iteration: 2000 / 10000 20% Iteration: 3000 / 10000 30% Iteration: 4000 / 10000 40% Iteration: 5000 / 10000 50% Iteration: 5001 / 10000 50% Iteration: 6000 / 10000 60% Iteration: 7000 / 10000 70% Iteration: 8000 / 10000 80% Iteration: 9000 / 10000 90% Iteration: 10000 / 10000 100%

Elapsed Time: 0.088028 seconds (Warm-up)

0.091544 seconds (Sampling)

0.179572 seconds (Total)

print(samples) Inference for Stan model: Rate_5_model. 1 chains, each with iter=10000; warmup=5000; thin=1; post-warmup draws per chain=5000, total post-warmup draws=5000.

         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat

theta 0.50 0.00 0.10 0.31 0.43 0.50 0.57 0.70 1632 1 postpredk1 5.02 0.03 1.88 1.00 4.00 5.00 6.00 9.00 3173 1 postpredk2 4.99 0.03 1.86 1.00 4.00 5.00 6.00 9.00 3026 1 lp__ -15.74 0.01 0.67 -17.65 -15.90 -15.47 -15.30 -15.25 2348 1

Samples were drawn using NUTS(diag_e) at Sat Dec 20 18:23:40 2014. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

theta <- extract(samples)$theta postpredk1 <- extract(samples)$postpredk1 postpredk2 <- extract(samples)$postpredk2

Two-panel plot.

layout(matrix(c(1,2),1,2)) layout.show(2)

First, a histogram for theta.

par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,

  • font.lab = 2, cex.axis = 1.3, bty = "n", las=1) Nbreaks <- 80 y <- hist(theta, Nbreaks, plot=F) plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
  • xlim=c(0,1), ylim=c(0,10), xlab="Theta", ylab="Density")

    let's plot a density estimate over this:

    lines(density(theta), col="red", lwd=2)

    Second plot, the data space (predictives)

plot(k1,k2,type="p", pch=4, cex=2, lwd=2, xlab="Success Count 1",

  • ylab="Success Count 2", xlim=c(-1, n1+1), ylim=c(-1,n2+1))
    Error in plot(k1, k2, type = "p", pch = 4, cex = 2, lwd = 2, xlab = "Success Count 1", : error in evaluating the argument 'x' in selecting a method for function 'plot': Error: object 'k1' not found nsamples <- length(theta) sc <- 10 for (i in 0:n1)
  • {
  • for (j in 0:n2)
  • {
  • match.preds <- sum(postpredk1==i & postpredk2==j)/nsamples
  • if (match.preds > 0)
  • {
  • points(i,j, pch=0, cex=sc*sqrt(match.preds))
  • }
  • }
  • } Error: object 'n1' not found
martin-smira commented 9 years ago

It is fixed now. Thanks again for reporting the issue.

drewyallop commented 9 years ago

Thank you!

Drew