Bnormal {bmstdr} | R Documentation |
N(theta, sigma2): Using different methods.
Description
N(theta, sigma2): Using different methods.
Usage
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:
|
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. |
Value
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)
}