######################################################################
# STAT 7630 Bayesian Statistics 
# Peng Zeng @ Auburn University
# 08-28-2025
######################################################################

######################################################################
# Markov chain
######################################################################

P = matrix(0, nrow = 6, ncol = 6)
P[c(1, 7, 8, 15, 22, 29, 30, 36)] = 0.5
P[c(2, 9, 14, 16, 21, 23, 28, 35)] = 0.25

size = 50000; status = numeric(size);
status[1] = 3
for(i in 2:size)
{
    status[i] = sample(1:6, size = 1, prob = P[status[i - 1], ])
}

table(status) / size

######################################################################
# mixture of normals  
######################################################################

mix.iid = function(n, prob, mu, sigma)
{
    k = length(mu)
    x = sample(1:k, size = n, prob = prob, replace = TRUE)
    rnorm(n, mu[x], sigma[x])
}

y1 = mix.iid(1000, c(0.45, 0.1, 0.45), c(-3, 0, 3), sqrt(c(1/3, 1/3, 1/3)))

par(mar = c(3, 3, 0.5, 0.5))
hist(y1, breaks = 30, prob = TRUE, xlab = "", main = "")
xgrid = seq(-6, 6, length = 201)
den = ( 0.45 * dnorm(xgrid, -3, sqrt(1/3)) 
      + 0.1  * dnorm(xgrid,  0, sqrt(1/3)) 
      + 0.45 * dnorm(xgrid,  3, sqrt(1/3)))
lines(xgrid, den, col = "red")


mix.Gibbs = function(n, prob, mu, sigma)
{
    k = length(mu)
    x = sample(1:k, size = 1, prob = prob) 
    y = rnorm(1, mu[x], sigma[x])
      
    sample.y = numeric(n); 
    for(i in 1:n)  {
        prob2 = prob * dnorm(y, mu, sigma)
        x = sample(1:k, size = 1, prob = prob2)
        y = rnorm(1, mu[x], sigma[x])
        sample.y[i] = y; 
    }
    sample.y
}

y2 = mix.Gibbs(1000, c(0.45, 0.1, 0.45), c(-3, 0, 3), sqrt(c(1/3, 1/3, 1/3)))

par(mar = c(3, 3, 0.5, 0.5))
hist(y2, breaks = 30, prob = TRUE, xlab = "", main = "")
xgrid = seq(-6, 6, length = 201)
den = ( 0.45 * dnorm(xgrid, -3, sqrt(1/3)) 
      + 0.1  * dnorm(xgrid,  0, sqrt(1/3)) 
      + 0.45 * dnorm(xgrid,  3, sqrt(1/3)))
lines(xgrid, den, col = "red")

par(mar = c(3, 3, 0.5, 0.5))
plot(y2, type = "l", xlab = "", ylab = "")

y3 = mix.Gibbs(100000, c(0.45, 0.1, 0.45), c(-3, 0, 3), sqrt(c(1/3, 1/3, 1/3)))

par(mar = c(3, 3, 0.5, 0.5))
hist(y3, breaks = 30, prob = TRUE, xlab = "", main = "")
xgrid = seq(-6, 6, length = 201)
den = ( 0.45 * dnorm(xgrid, -3, sqrt(1/3)) 
      + 0.1  * dnorm(xgrid,  0, sqrt(1/3)) 
      + 0.45 * dnorm(xgrid,  3, sqrt(1/3)))
lines(xgrid, den, col = "red")

######################################################################
# sampling from t-distribution 
# with normal distribution as the proposal distribution 
######################################################################

rt.RW = function(n, tdf, x0, sigma)
{
    x = numeric(n)
    x[1] = x0
    for(i in 2:n)
    {
        y = rnorm(1, x[i-1], sigma)
        logratio = (dt(y,      df = tdf, log = TRUE)
                  - dt(x[i-1], df = tdf, log = TRUE))
        x[i] = ifelse(log(runif(1)) <= logratio, y, x[i-1])
    }
    x
}

x1 = rt.RW(1000, 4, 25, 0.05)
x2 = rt.RW(1000, 4, 25, 0.5)
x3 = rt.RW(1000, 4, 25, 2)
x4 = rt.RW(1000, 4, 25, 16)

par(mfrow = c(2, 2), mar = c(3, 3, 0.5, 0.5))
plot(x1, type = "l")
plot(x2, type = "l")
plot(x3, type = "l")
plot(x4, type = "l")

acf(x1)
acf(x2)
acf(x3)
acf(x4)

mean(diff(x1) != 0)
mean(diff(x2) != 0)
mean(diff(x3) != 0)
mean(diff(x4) != 0)

######################################################################
# sampling from Rayleigh distribution
# with chi-squared distribution as the proposal distribution
######################################################################

logdRay = function(x, sigma)
{  log(x) - 2 * log(sigma) - x * x / (2 * sigma^2)  }

rRay = function(n, sigma = 1)
{
    x = numeric(n)
    x[1] = sigma 
    for(i in 2:n)
    {
        y = rchisq(1, df = x[i-1])
        logratio = (logdRay(y, sigma) - logdRay(x[i-1], sigma)
                    + dchisq(x[i-1], df = y, log = TRUE)
                    - dchisq(y, df = x[i-1], log = TRUE))
        x[i] = ifelse(log(runif(1)) <= logratio, y, x[i-1])
    }
    x
}

x = rRay(500, sigma = 4);

par(mfrow = c(1, 3), mar = c(3, 3, 0.5, 0.5))
plot(x, type = "l")
acf(x)
hist(x, breaks = 50, prob = TRUE, main = "")
xseq = seq(0, max(x), length = 201);
lines(xseq, exp(logdRay(xseq, 4)), col = "red")

######################################################################
# height of men 
######################################################################

d = list(age = c(-Inf, 66, 68, 70, 72, 74, Inf),
        freq = c(14, 30, 49, 70, 33, 15))

logden = function(theta, mylist)
{
    mu = theta[1]; 
    sigma = exp(theta[2]);
    prob = pnorm(mylist$age, mu, sigma)
    sum(mylist$freq * log(diff(prob)));
}

x = rep(seq(65, 75, by = 2), c(14, 30, 49, 70, 33, 15))
c(mean(x), log(sd(x)))

fit = optim(c(70.17, 0.95), logden, hessian = TRUE,
            mylist = d, control = list(fnscale = -1));
V = solve(-fit$hessian)

size = 5000
theta = c(70, 1); 
thmat = matrix(0, nrow = size, ncol = 2);
count = 0; 
for(i in 1:size)
{
    theta.new = theta + 2 * MASS::mvrnorm(1, c(0, 0), V); 
    logratio = (logden(theta.new, d) - logden(theta, d));
    if(log(runif(1)) < logratio) {
        count = count + 1
        theta = theta.new;
    } 
    
    thmat[i, ] = theta; 
}

count / size 
mean(thmat[-c(1:2000), 1])
mean(exp(thmat[-c(1:2000), 2]))
quantile(thmat[-c(1:2000), 1], prob = c(0.025, 0.5, 0.975))
quantile(exp(thmat[-c(1:2000), 2]), prob = c(0.025, 0.5, 0.975))

par(mfrow = c(2, 2), mar = c(3, 3, 0.5, 0.5))
plot(thmat[-(1:2000), 1], type = "l", xlab = "")
plot(thmat[-(1:2000), 2], type = "l", xlab = "")
acf(thmat[-(1:2000), 1], main = "", xlab = "")
acf(thmat[-(1:2000), 2], main = "", xlab = "")

######################################################################
# height of men - stan 
######################################################################

library(rstan)

data_list = list(K = 6, 
       cuts = c(66, 68, 70, 72, 74),
       counts = c(14, 30, 49, 70, 33, 15))

fit = stan(
    file = "grouped-normal.stan",    # Stan program
    data = data_list,                # named list of data
    init = list(list(mu = 70, sigma = 1), 
                list(mu = 65, sigma = 1.5), 
                list(mu = 75, sigma = 0.8), 
                list(mu = 70, sigma = 2)),
    chains = 4,             # # of Markov chains
    warmup = 2000,          # # of warmup iterations per chain
    iter = 5000,            # # of iterations per chain
    cores = 1,              # number of cores per chain
    refresh = 0             # no progress shown
)

print(fit)
traceplot(fit, pars = c("mu", "sigma"))

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