Bnormal function

N(theta, sigma2): Using different methods.

Bnormal(
  package = "exact",
  y = ydata,
  mu0 = mean(y),
  kprior = 0,
  prior.M = 1e-04,
  prior.sigma2 = c(0, 0),
  N = 2000,
  burn.in = 1000,
  rseed = 44
)

Arguments

• package: Which package (or method) to use. Possibilities are:

  • "exact": Use exact theoretical calculation.
  • "RGibbs": Use Gibbs sampler using R code.
  • "stan": Use HMC by implementing in Stan.
  • "inla": Use the INLA package.

• y: A vector of data values. Default is 28 ydata values from the package bmstdr

• mu0: The value of the prior mean if kprior=0. Default is the data mean.

• kprior: A scalar providing how many data standard deviation the prior mean is from the data mean. Default value is 0.

• prior.M: Prior sample size, defaults to 10^(-4).#'

• prior.sigma2: Shape and scale parameter value for the gamma prior on 1/sigma^2, the precision.

• N: is the number of Gibbs sampling iterations

• burn.in: is the number of initial iterations to discard before making inference.

• rseed: is the random number seed defaults to 44.

Returns

A list containing the exact posterior means and variances of theta and sigma2

Examples

# Find the posterior mean and variance  using `exact' methods - no sampling 
# or approximation  
Bnormal(kprior = 1, prior.M = 1, prior.sigma2 = c(2, 1))
# Use default non-informative prior
Bnormal(mu0 = 0) 
# Start creating table
y <-  ydata
mu0 <-  mean(y)
kprior <-  1
prior.M  <-  1
prior.sigma2 <- c(2, 1)
N  <-  10000
eresults <- Bnormal(package = "exact", y = y, mu0 = mu0, kprior = kprior,
    prior.M = prior.M, prior.sigma2 = prior.sigma2)
eresults 
# Run Gibbs sampling
samps <- Bnormal(package = "RGibbs", y = y, mu0 = mu0, kprior = kprior,
    prior.M = prior.M, prior.sigma2 = prior.sigma2, N = N)
gres <- list(mean_theta = mean(samps[, 1]), var_theta = var(samps[, 1]),
    mean_sigma2 = mean(samps[, 2]), var_sigma2 = var(samps[, 2]))
glow <- list(theta_low = quantile(samps[, 1], probs = 0.025), var_theta = NA,
    sigma2_low = quantile(samps[, 2], probs = 0.025), var_sigma2 = NA)
gupp <- list(theta_low = quantile(samps[, 1], probs = 0.975), var_theta = NA,
    sigma2_low = quantile(samps[, 2], probs = 0.975), var_sigma2 = NA)
a <- rbind(unlist(eresults), unlist(gres), unlist(glow), unlist(gupp))
cvsamp <- sqrt(samps[, 2])/samps[, 1]
cv <- c(NA, mean(cvsamp), quantile(cvsamp, probs = c(0.025, 0.975)))
u <- data.frame(a, cv)
rownames(u) <- c("Exact", "Estimate", "2.5%", "97.5%")
print(u)
# End create table 
##
# Compute using the model fitted by Stan 
u <- Bnormal(package = "stan", kprior = 1, prior.M = 1, prior.sigma = c(2, 1),
    N = 2000, burn.in = 1000)
print(u)
###
# Compute using INLA 
if (require(INLA)) {
    u <- Bnormal(package = "inla", kprior = 1, prior.M = 1, prior.sigma = c(2,
        1), N = 1000)
    print(u)
}