######################################################################
# STAT 7630 Bayesian Statistics 
# Peng Zeng @ Auburn University
# 09-04-2025
######################################################################

######################################################################
# soccer goals 
######################################################################

n = 35; s = 57;

logpost1 = function(loglda, hyper.pars)
{
    lambda = exp(loglda); 
    (dgamma(lambda, s + 1, n, log = TRUE) 
        + dgamma(lambda, hyper.pars[1], hyper.pars[2], log = TRUE) 
        + log(lambda)); 
}

fit1 = optim(log(s/n), logpost1, hyper.pars = c(4.57, 1.43), 
        lower = -log(s), upper = log(s), 
        method = "Brent", hessian = TRUE, control = list(fnscale = -1));
x1 = fit1$par;  
mar1 = sqrt(2*pi) / sqrt(-fit1$hessian) * exp(fit1$value);

logpost2 = function(loglda, hyper.pars)
{
    lambda = exp(loglda); 
    (dgamma(lambda, s + 1, n, log = TRUE) 
        + dnorm(loglda, hyper.pars[1], hyper.pars[2], log = TRUE)); 
}

fit2 = optim(log(s/n), logpost2, hyper.pars = c(1, 0.5), 
        lower = -log(s), upper = log(s), 
        method = "Brent", hessian = TRUE, control = list(fnscale = -1));
x2 = fit2$par; 
mar2 = sqrt(2*pi) / sqrt(-fit2$hessian) * exp(fit2$value);

fit3 = optim(log(s/n), logpost2, hyper.pars = c(2, 0.5), 
        lower = -log(s), upper = log(s), 
        method = "Brent", hessian = TRUE, control = list(fnscale = -1));
x3 = fit3$par; 
mar3 = sqrt(2*pi) / sqrt(-fit3$hessian) * exp(fit3$value);

fit4 = optim(log(s/n), logpost2, hyper.pars = c(1, 2), 
        lower = -log(s), upper = log(s), 
        method = "Brent", hessian = TRUE, control = list(fnscale = -1));
x4 = fit4$par; 
mar4 = sqrt(2*pi) / sqrt(-fit4$hessian) * exp(fit4$value);

######################################################################
# THE END
######################################################################