rmcelreath

\newwrite\tempfile

\usepackage{verbatim}
\makeatletter
\newwrite\Code@out % temp file for writing out and reading back in for display

\newcommand\VerbSaver{\obeylines\expandafter\VerbSaverArg\noexpand}

\newcommand\VerbSaverArg[1][code.txt]{%
    \gdef\FName{xcodetempx.txt}%

rmcelreath

######################
# Replication and reanalysis for:
#   van der Lee and Ellemers 10.1073/pnas.1510159112

# load data from supplemental Table S1
# found at http://www.pnas.org/content/early/2015/09/16/1510159112/suppl/DCSupplemental

d <- read.csv("http://xcelab.net/VDL_data.csv")
d$rejects <- d$applications - d$awards

rmcelreath

# functions for plotting garden of forking data plots

library(rethinking)

polar2screen <- function( dist, origin, theta ) {
  ## takes dist, angle and origin and returns x and y of destination point
  vx <- cos(theta) * dist;
  vy <- sin(theta) * dist;
  c( origin[1]+vx , origin[2]+vy );
}

rmcelreath

# "impute" missing binary predictor
# really just marginalizes over missingness
# imputed values produced in generated quantities

N <- 1000 # number of cases
N_miss <- 100 # number missing values
x_baserate <- 0.25 # prob x==1 in total sample
a <- 0 # intercept in y ~ N( a+b*x , 1 )
b <- 1 # slope in y ~ N( a+b*x , 1 )

rmcelreath

library(rstan)

# simple model with a bug
model_code <- "
parameters{
    real a;
}
model{
    a ~ normal(0,1);
}

rmcelreath

# berkson's paradox journal example

library(rethinking)

n <- 500
rho <- 0.25
y <- rmvnorm2( n , Mu=c(0,0) , sigma=c(1,1) , Rho=matrix( c(1,rho,rho,1) , 2 , 2 ) )
b <- 1
score <- y[,1] + b*y[,2]
theshold <- quantile( score , 0.9 )

rmcelreath

# script for examples in Bayes@Lund2017 presentation

# joint model example
notes_max <- 10
rate_max <- 5
pm <- matrix( NA , nrow=rate_max+1 , ncol=notes_max+1 )
for ( i in 1:(rate_max+1) )
    for ( j in 1:(notes_max+1) )
        pm[i,j] <- dpois( j-1 , lambda=i ) * (1/(rate_max+1))

rmcelreath

# chapter 2 code

library(rethinking)

col1 <- "black"
col2 <- col.alpha( "black" , 0.7 )

#####
# introducing the golem
# show posterior as data comes in

rmcelreath

p <- list()
p$A <- c(0,0,10,0,0)
p$B <- c(0,1,8,1,0)
p$C <- c(0,2,6,2,0)
p$D <- c(1,2,4,2,1)
p$E <- c(2,2,2,2,2)

p_norm <- lapply( p , function(q) q/sum(q))
( H <- sapply( p_norm , function(q) -sum(ifelse(q==0,0,q*log(q))) ) )

rmcelreath

library(rethinking)
p_grid <- seq( from=0 , to=1 , length.out=100 )
prior <- rep( 1 , 100 )
likelihood <- dbinom( 6 , size=9 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
plot( posterior ~ p_grid , type="l" )
shade( posterior ~ p_grid , lim=c(0,0.5) , col=rangi2 )