# Testing software instalations
# author: Peter McMahan
# date: 2024-12-26

# load packages
library(rethinking)
library(brms)
library(lme4)

# draw 10,000 random normal samples and look at them using dens()
obs <- rnorm(10e3,mean = 2,sd = 1) 
dens(obs)
dens(obs,norm.comp = TRUE)



###
# Testing rethinking functions
###

# test quap() function:
# make a simple model
f <- alist(
  y ~ dnorm( mu , sigma ),
  mu ~ dnorm( 0 , 20 ),
  sigma ~ dcauchy( 0 , 3 )
)
# fit it using maximum a posteriori / quadratic approximation
fit_quap <- quap( 
  f , 
  data = list(y=obs) ,
  start=list(mu=0,sigma=1)
)
summary(fit_quap)


# test the ulam() function
# (this should take a while to run, and will
# output some text while it does so. Be patient!)
# This fits the same model as above, using MCMC instead of MAP
fit_stan <- ulam( 
  f , 
  data = list(y=obs)
)
summary(fit_stan)



###
# Testing brms
###
fit_brms <- brm(y ~ 1, data=list(y=obs))


###
# Testing lme4
###
fit_lme4 <- lmer(y ~ 1 + (1 | c), data=data.frame(y=obs, c=c(1,2)))



###
# compare the outputs, for fun
# (the intercept estimation should be nearly
# identical for all models)
###
print(precis(fit_quap))
print(precis(fit_stan))
print(summary(fit_brms))
print(summary(fit_lme4))