search

stan-dev / example-models

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

Survey_Stan.R #14

Closed drewyallop closed 9 years ago

drewyallop commented 9 years ago

Does not plot. Console ...

clears workspace:

rm(list=ls())

library(rstan) Loading required package: Rcpp Loading required package: inline

Attaching package: ‘inline’

The following object is masked from ‘package:Rcpp’:

registerPlugin

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

to be passed on to Stan

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

myinits <- list(

  • list(theta=.5))

    parameters to be monitored:

parameters <- c("theta", "n")

The following command calls Stan with specific options.

For a detailed description type "?stan".

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

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

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

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

  • )

TRANSLATING MODEL 'Survey_model' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'Survey_model' NOW. In file included from file29b37e27a154.cpp:466: 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 file29b37e27a154.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 file29b37e27a154.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 'Survey_model' NOW (CHAIN 1).

Iteration: 1 / 4000 0% Iteration: 400 / 4000 10% Iteration: 800 / 4000 20% Iteration: 1200 / 4000 30% Iteration: 1600 / 4000 40% Iteration: 2000 / 4000 50% Iteration: 2001 / 4000 50% Iteration: 2400 / 4000 60% Iteration: 2800 / 4000 70% Iteration: 3200 / 4000 80% Iteration: 3600 / 4000 90% Iteration: 4000 / 4000 100%

Elapsed Time: 13.5233 seconds (Warm-up)

15.1188 seconds (Sampling)

28.6422 seconds (Total)

Now the values for the monitored parameters are in the "samples" object,

ready for inspection.

theta <- extract(samples)$theta n <- extract(samples)$n

First calculate MLE:

cc <- -Inf ind <- 0

for (i in 1:length(n)) {

  • logL <- 0
  • for(j in 1:m) {
  • logL <- logL+lgamma(n[i]+1)-lgamma(k[j]+1)-lgamma(n[i]-k[j]+1)
  • logL <- logL+k[j]_log(theta[i])+(n[i]-k[j])_log(1-theta[i])
  • }
  • if (logL>cc) {
  • ind <- i
  • cc <- logL
  • }
  • } Error: object 'm' not found

    end MLE

    Plots

layout(matrix(c(2,0,1,3),2,2,byrow=T), width=c(2/3, 1/3), heights=c(1/3, 2/3)) xhist <- hist(n, plot=F) yhist <- hist(theta, plot=F) top <- max(c(xhist$counts, yhist$counts)) xrange <- c(0, nmax) Error: object 'nmax' not found yrange <- c(0, 1)

par(mar=c(5,5,1,1)) plot(n,theta,xlim=xrange, ylim=yrange,ylab="", xlab="") Error in plot.default(n, theta, xlim = xrange, ylim = yrange, ylab = "", : object 'xrange' not found axis(1) Error in axis(1) : plot.new has not been called yet mtext("Number of Surveys", side=1,line=2.25, cex=1.2) Error in mtext("Number of Surveys", side = 1, line = 2.25, cex = 1.2) : plot.new has not been called yet axis(2, cex=1.2) Error in axis(2, cex = 1.2) : plot.new has not been called yet las=0 mtext("Rate of Return", side=2 ,line=2.25, cex=1.2) Error in mtext("Rate of Return", side = 2, line = 2.25, cex = 1.2) : plot.new has not been called yet las=1 points(mean(n),mean(theta), col="red", lwd=3, pch=4) #expectation Error in plot.xy(xy.coords(x, y), type = type, ...) : plot.new has not been called yet points(n[ind],theta[ind], col="green", lwd=3, pch=10) #Maximum Likelihood Error in plot.xy(xy.coords(x, y), type = type, ...) : plot.new has not been called yet

par(mar=c(0,4,1,1)) barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0,col="lightblue")

par(mar=c(4,0,1,3)) barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE,

  • col="lightblue")
martin-smira commented 9 years ago

Thanks again. Sorry about these issues, in the further models you won't find so many issues. Here I separated the model code, data, and the R script in separate files and forgot to link the data correctly in the R script. It is fixed now!

drewyallop commented 9 years ago

Thank you. Drew

Sent from my iPad

On Dec 22, 2014, at 16:54, neekro notifications@github.com wrote:

Thanks again. Sorry about these issues, in the further models you won't find so many issues. Here I separated the model code, data, and the R script in separate files and forgot to link the data correctly in the R script. It is fixed now!

— Reply to this email directly or view it on GitHub.