rmcelreath

library(rethinking)

HMC2 <- function (U, grad_U, epsilon, L, current_q , ... ) {
  q = current_q
  p = rnorm(length(q),0,1)  # independent standard normal variates
  current_p = p
  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U( q , ... ) / 2
  # Alternate full steps for position and momentum
  qtraj <- matrix(NA,nrow=L+1,ncol=length(q))

rmcelreath

# script to illustrate collider bias described in Breen 2018 doi:10.1093/esr/jcy037

library(rethinking)

N <- 200  # number of grandparent-parent-child triads
b_gp <- 1 # direct effect of G on P
b_gc <- 0 # direct effect of G on C
b_pc <- 1 # direct effect of P on C
b_U <- 2  # direct effect of U on P and C

rmcelreath

// simple reinforcement model
data{
    int N;          // total number of observations
    int N_trials;   // num trials per individual
    int N_individuals;
    int trial[N];   // trial number in individual
    int y[N];       // choices
    real P[N,2];       // payoffs to each option
    int id[N];      // individual ID
}

rmcelreath

// contest model
data{
    int N_choices;
    int N_items;
    int N_raters;
    int y[N_choices];
    int item1[N_choices];
    int item2[N_choices];
    int rater[N_choices];
    real x[N_items];      // item feature

rmcelreath

###################
# example workflow
# 1 build model
# 2 simulate data
# 3 analyze data

################
# 1 build model

# example context: tradeoff between egg number and egg survival

rmcelreath

## code to reproduce Figure 2.5 on page 30 of Statistical Rethinking

library(rethinking)

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

# show posterior as data comes in
# 1 indicates 'water'; 0 indicates 'land'
d <- c(1,0,1,1,1,0,1,0,1) # length 9

rmcelreath

# sim confounding of age --> happiness by controlling for marital status, when in truth age + happiness --> marriage
# assume:
# (1) no relationship between happiness and age
# (2) 5 happiest people each year get married (after age 20)

N_years <- 200
N_born_per_year <- 40
N_marry_per_year <- 5
mortality <- 0.05
max_age <- 100

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 )

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

# chapter 2 code

library(rethinking)

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

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