Return to this later
rm(list=ls())
set.seed(874)
startTime<-proc.time()
library(knitr)
options(width=120)
opts_chunk$set(comment = "", warning = FALSE, message = FALSE,
echo = TRUE, tidy = FALSE, size="tiny", cache=FALSE,
progress=TRUE,
cache.path = 'program_Cache/',
fig.path='figure/')
knitr::knit_hooks$set(inline = function(x) {
knitr:::format_sci(x, 'md')
}) where<-"home" #this is used in the sourced program
path <- ""
work<- paste("X:/", path, sep = "")
nonwork<- paste("~/X/", path, sep = "")
if (where=="home") {wd<- nonwork} else {wd<-work}
path2 <- ""
work2<- paste("X:/", path2, sep = "")
nonwork2<- paste("~/X/", path2, sep = "")
if (where=="home") {wd2<- nonwork2} else {wd2<-work2}
work3<- paste("X:/FUNCTIONS/R", sep = "")
nonwork3<- paste("~/X/FUNCTIONS/R", sep = "")
if (where=="home") {wd3<- nonwork3} else {wd3<-work3}
setwd(wd)
opts_knit$set(root.dir = wd) ##THIS SETS YOUR WORKING DIRECTORY list.of.packages <- c("lme4","rethinking", "rstanarm")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
sapply(X = list.of.packages, require, character.only = TRUE)
p3 <- function(x) {formatC(x, format="f", digits=3)}
p4 <- function(x) {formatC(x, format="f", digits=4)}
p2 <- function(x) {formatC(x, format="f", digits=2)}
p1 <- function(x) {formatC(x, format="f", digits=0)}
# p1 <- function(x) {print(formatC(x, format="f", digits=1),quote=FALSE)}
# p2 <- function(x) {print(formatC(x, format="f", digits=2),quote=FALSE)}
# p3 <- function(x) {print(formatC(x, format="f", digits=3),quote=FALSE)}
# p4 <- function(x) {print(formatC(x, format="f", digits=4),quote=FALSE)}
#perhaps help colour plot text based on loop count
is.even <- function(x){ x %% 2 == 0 } data(eggs, package="faraway")
eggs$labtech <- factor(paste0(eggs$Lab,eggs$Technician))
eggs$labtechsamp <- factor(paste0(eggs$Lab,eggs$Technician,eggs$Sample))
levind1 <- as.numeric(eggs$Lab)
levind2 <- as.numeric(eggs$labtech)
levind3 <- as.numeric(eggs$labtechsamp)
sdscal <- sd(eggs$Fat)
eggdat <- list(Nobs=nrow(eggs),
Nlev1=max(levind1),
Nlev2=max(levind2),
Nlev3=max(levind3),
y=eggs$Fat,
levind1=levind1,
levind2=levind2,
levind3=levind3,
sdscal=sdscal)
nested.stan <-"
data {
int<lower=0> Nobs;
int<lower=0> Nlev1;
int<lower=0> Nlev2;
int<lower=0> Nlev3;
vector[Nobs] y;
int<lower=1,upper=Nlev1> levind1[Nobs];
int<lower=1,upper=Nlev2> levind2[Nobs];
int<lower=1,upper=Nlev3> levind3[Nobs];
real<lower=0> sdscal;
}
parameters {
real mu;
real<lower=0> sigmalev1;
real<lower=0> sigmalev2;
real<lower=0> sigmalev3;
real<lower=0> sigmaeps;
vector[Nlev1] eta1;
vector[Nlev2] eta2;
vector[Nlev3] eta3;
}
transformed parameters {
vector[Nlev1] ran1;
vector[Nlev2] ran2;
vector[Nlev3] ran3;
vector[Nobs] yhat;
ran1 = sigmalev1 * eta1;
ran2 = sigmalev2 * eta2;
ran3 = sigmalev3 * eta3;
for (i in 1:Nobs)
yhat[i] <- mu+ran1[levind1[i]]+ran2[levind2[i]]+ran3[levind3[i]];
}
model {
eta1 ~ normal(0, 1);
eta2 ~ normal(0, 1);
eta3 ~ normal(0, 1);
sigmalev1 ~ cauchy(0, 2.5*sdscal);
sigmalev2 ~ cauchy(0, 2.5*sdscal);
sigmalev3 ~ cauchy(0, 2.5*sdscal);
sigmaeps ~ cauchy(0, 2.5*sdscal);
y ~ normal(yhat, sigmaeps);
}"
rt <- stanc( model_code="nested.stan", model_name="Nested")
sm <- stan_model(stanc_ret = rt, verbose=FALSE)C:/Users/User/Documents/R/win-library/3.2/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function] system.time(fit <- sampling(sm, data=eggdat))
SAMPLING FOR MODEL 'Nested' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.734 seconds (Warm-up)
0.752 seconds (Sampling)
2.486 seconds (Total)
[1] "The following numerical problems occured the indicated number of times after warmup on chain 1"
count
Exception thrown at line 46: normal_log: Location parameter[1] is 1.#INF, but must be finite! 1
[1] "When a numerical problem occurs, the Metropolis proposal gets rejected."
[1] "However, by design Metropolis proposals sometimes get rejected even when there are no numerical problems."
[1] "Thus, if the number in the 'count' column is small, do not ask about this message on stan-users."
SAMPLING FOR MODEL 'Nested' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.542 seconds (Warm-up)
0.822 seconds (Sampling)
2.364 seconds (Total)
SAMPLING FOR MODEL 'Nested' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.508 seconds (Warm-up)
0.868 seconds (Sampling)
2.376 seconds (Total)
SAMPLING FOR MODEL 'Nested' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 1.485 seconds (Warm-up)
0.755 seconds (Sampling)
2.24 seconds (Total) user system elapsed
9.56 0.02 9.61 print(fit,pars=c("mu","sigmalev1","sigmalev2","sigmalev3","sigmaeps"))Inference for Stan model: Nested.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu 0.39 0 0.06 0.25 0.35 0.39 0.42 0.50 481 1.00
sigmalev1 0.09 0 0.06 0.01 0.05 0.08 0.12 0.26 644 1.01
sigmalev2 0.09 0 0.04 0.01 0.06 0.09 0.12 0.19 861 1.00
sigmalev3 0.06 0 0.03 0.00 0.04 0.06 0.08 0.12 515 1.01
sigmaeps 0.09 0 0.01 0.07 0.08 0.09 0.10 0.12 1359 1.00
Samples were drawn using NUTS(diag_e) at Thu Jul 21 22:31:58 2016.
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).
The estimated Bayesian Fraction of Missing Information is a measure of
the efficiency of the sampler with values close to 1 being ideal.
For each chain, these estimates are
0.6 0.7 0.7 0.7R version 3.2.2 (2015-08-14)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8 x64 (build 9200)
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] rstanarm_2.10.1 Rcpp_0.12.4 rethinking_1.58 rstan_2.10.1 StanHeaders_2.10.0-2
[6] ggplot2_2.1.0 lme4_1.1-12 Matrix_1.2-2 knitr_1.12.3
loaded via a namespace (and not attached):
[1] gtools_3.5.0 zoo_1.7-13 shinyjs_0.6 reshape2_1.4.1 splines_3.2.2 lattice_0.20-33
[7] colorspace_1.2-6 miniUI_0.1.1 htmltools_0.3.5 stats4_3.2.2 loo_0.1.6 yaml_2.1.13
[13] base64enc_0.1-3 nloptr_1.0.4 matrixStats_0.50.2 plyr_1.8.3 stringr_1.0.0 munsell_0.4.3
[19] gtable_0.2.0 htmlwidgets_0.6 mvtnorm_1.0-5 codetools_0.2-14 threejs_0.2.2 coda_0.18-1
[25] evaluate_0.9 inline_0.3.14 httpuv_1.3.3 markdown_0.7.7 xts_0.9-7 xtable_1.8-2
[31] scales_0.4.0 DT_0.1 shinystan_2.2.0 formatR_1.3 jsonlite_0.9.19 mime_0.4
[37] gridExtra_2.2.1 digest_0.6.9 stringi_1.0-1 shiny_0.13.2 grid_3.2.2 tools_3.2.2
[43] magrittr_1.5 shinythemes_1.0.1 MASS_7.3-45 rsconnect_0.4.3 dygraphs_0.9 minqa_1.2.4
[49] rmarkdown_0.9.6 R6_2.1.2 nlme_3.1-127 [1] "~/X/"This took 3.79 seconds to execute.