rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Greg Francis
# PSY 626
# 31 October 2025


# Bayes Factor 

library(BayesFactor)

scores<- c(12, 10, 9, 13, 13, 8, 11, 7, 14, 11, 15, 17, 11, 12)

# Subtract the value in the null hypothesis (create difference scores)

diffScores = scores - 10

bf = ttestBF(x = diffScores, nullInterval = c(-Inf, 0))

bf

print(bf[2]/bf[1])

# Stan/Ulam with rethinking 

# set up data frame

data <- data.frame(Participant=seq(1,14),  Score =scores )		

library(rethinking)

# Model comparison
# Null model
NullModel <- ulam(
	alist( Score ~ dnorm(10, sigma), 
	sigma ~ dnorm(0, 10)
	), data= data, constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 

# Alt model
AltModel <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu ~ dnorm(0, 10),
	sigma ~ dnorm(0, 10)
	), data= data,  constraints=list(mu="lower=10,upper=100", sigma="lower=0"), log_lik=TRUE, chains=4  ) 

compare(NullModel, AltModel)


# Analyzing posterior

# Full model
FullModel <- ulam(
	alist( Score ~ dnorm(mu, sigma), 
	mu ~ dnorm(0, 10),
	sigma ~ dnorm(0, 10)
	), data= data,  constraints=list(sigma="lower=0"), log_lik=TRUE, chains=4  ) 
	
#plot posterior of a
numSamples=2000
post<-extract.samples(FullModel, n= numSamples)

dev.new()
dens(post$mu, col=rangi2, lwd=2, xlab="posterior for mu")

# Probability that mu > 10 ?
length(post$mu[post$mu >10])/length(post$mu)