Use naive bayes to calculate distribution mean + credible interval

iris_naive_bayes.R

# http://math.utoledo.edu/~mleite/math-EES-seminar/week1Baysean.pdf

data(iris)

# Estimation ####

y     <- iris$Sepal.Length
y_bar <- mean(y)
y_sd  <- sd(y)
n     <- length(y)

# Confidence Interval ####

ci_mu    <- y_bar + y_sd/sqrt(n) * qt(c(0.025, 0.975), n - 1)
ci_sigma <- y_sd * sqrt((n - 1)) / qchisq(c(0.025, 0.975), n - 1)

# Naive Bayes ####

n.sims <- 5000
chi2 <- rchisq(n.sims, df = n - 1)  
sigma <- y_sd * sqrt((n - 1) / chi2) # likelihood * prior ?  
mu <- rnorm(n.sims, y_bar, sigma/sqrt(n))
y_tilde <- rnorm(n.sims, mu, sigma)

# Plotting ####
par(mfrow = c(1, 2))
hist(y)
abline(v = c(y_bar, ci_mu), col = "red")

hist(y_tilde)
abline(v = c(y_bar, mu[order(mu)][length(mu) * c(.025, 0.975 )]), col = "red") # mean +- 95% credible interval