rmcelreath

# medals per capita demonstration

d <- read.delim("medals_all.tsv")

library(rethinking)
# blank2()

d$Medals <- d$WeightedMedals
# d$Medals <- d$Gold + d$Silver + d$Bronze
d$Pop <- as.numeric(gsub(",", "", d$Population))

rmcelreath

library(rethinking)

set.seed(73)
N <- 500
U_sim <- rnorm( N )
Q_sim <- sample( 1:4 , size=N , replace=TRUE )
E_sim <- rnorm( N , U_sim + Q_sim )
W_sim <- rnorm( N , U_sim + 0*E_sim )

Rsize <- 2

rmcelreath

# prior - likelihood conflict

library(rethinking)

yobs <- 0

mtt <- ulam(
    alist(
        y ~ dstudent(2,mu,1),
        mu ~ dstudent(2,10,1)

rmcelreath

# distribution of null p-values in binomial glm

library(rethinking)

# t test
# pval has uniform distribution under null

f <- function(N=10) {
    y1 <- rnorm(N)
    y2 <- rnorm(N)

rmcelreath

# simulation of simple Kopernican model showing epicycle

library(rethinking)
library(animation)

my_circle <- function(x=0,y=0,r=1,angle=0,...) {
    a <- seq(angle, angle + 2 * pi, length = 360)
    lines( r*cos(a)+x , r*sin(a)+y , ... )
}
# plot(NULL,xlim=c(-1,1),ylim=c(-1,1)); my_circle(0,0,0.5,lty=2)

rmcelreath

library(rethinking)
library(animation)

ani.saveqz <- function (list , dpi=150) {   
    if (missing(list)) 
        list = animation:::.ani.env$.images
    lapply(1:length(list), function(x) {
        dev.hold()
        replayPlot( list[[x]] )
        ani.pause()

rmcelreath

# version that marginalizes out the missing data
flbi_plus <- ulam(
    alist(
        c(M,D) ~ multi_normal( c(mu,nu) , Rho , Sigma ),
        mu <- a1 + b*B1,
        nu <- a2 + b*B2 + m*M,
        c(a1,a2,b,m) ~ normal( 0 , 0.5 ),
        Rho ~ lkj_corr( 2 ),
        Sigma ~ exponential( 1 )
    ), data=dat , chains=4 , cores=4 , cmdstan=TRUE )

rmcelreath

# more exotic example - no instrument (B1 -> D), but have two measures of U
set.seed(1908)
N <- 200 # number of pairs
U <- rnorm(N,0,1) # simulate confound
V <- rnorm(N,U,1)
W <- rnorm(N,U,1)
# birth order and family sizes
B1 <- rbinom(N,size=1,prob=0.5) # 50% first borns
M <- rnorm( N , 2*B1 + U )
B2 <- rbinom(N,size=1,prob=0.5)

rmcelreath

# simulate interventions

# B1 -> D (should be b*m)
post <- extract.samples(flbi)
quantile( with(post,b*m) )

# now do the same with simulation

# B1 = 0
# B1 -> M

rmcelreath

# full-luxury bayesian inference
library(rethinking)
library(cmdstanr)
dat <- list(N=N,M=M,D=D,B1=B1,B2=B2)
set.seed(1908)
flbi <- ulam(
    alist(
        # mom model
            M ~ normal( mu , sigma ),
            mu <- a1 + b*B1 + k*U[i],