Binomial sampling with a general continuous prior
Evaluates and plots the posterior density for pi, the probability of a success in a Bernoulli trial, with binomial sampling and a general continuous prior on pi
binogcp(
x,
n,
density = c("uniform", "beta", "exp", "normal", "user"),
params = c(0, 1),
n.pi = 1000,
pi = NULL,
pi.prior = NULL,
...
)x |
the number of observed successes in the binomial experiment. |
n |
the number of trials in the binomial experiment. |
density |
may be one of "beta", "exp", "normal", "student", "uniform" or "user" |
params |
if density is one of the parameteric forms then then a vector of parameters must be supplied. beta: a, b exp: rate normal: mean, sd uniform: min, max |
n.pi |
the number of possible pi values in the prior |
pi |
a vector of possibilities for the probability of success in a single trial. This must be set if density = "user". |
pi.prior |
the associated prior probability mass. This must be set if density = "user". |
... |
additional arguments that are passed to |
A list will be returned with the following components:
likelihood |
the scaled likelihood function for pi given x and n |
posterior |
the posterior probability of pi given x and n |
pi |
the vector of possible pi values used in the prior |
pi.prior |
the associated probability mass for the values in pi |
## simplest call with 6 successes observed in 8 trials and a continuous
## uniform prior
binogcp(6, 8)
## 6 successes, 8 trials and a Beta(2, 2) prior
binogcp(6, 8,density = "beta", params = c(2, 2))
## 5 successes, 10 trials and a N(0.5, 0.25) prior
binogcp(5, 10, density = "normal", params = c(0.5, 0.25))
## 4 successes, 12 trials with a user specified triangular continuous prior
pi = seq(0, 1,by = 0.001)
pi.prior = rep(0, length(pi))
priorFun = createPrior(x = c(0, 0.5, 1), wt = c(0, 2, 0))
pi.prior = priorFun(pi)
results = binogcp(4, 12, "user", pi = pi, pi.prior = pi.prior)
## find the posterior CDF using the previous example and Simpson's rule
myCdf = cdf(results)
plot(myCdf, type = "l", xlab = expression(pi[0]),
ylab = expression(Pr(pi <= pi[0])))
## use the quantile function to find the 95% credible region.
qtls = quantile(results, probs = c(0.025, 0.975))
cat(paste("Approximate 95% credible interval : ["
, round(qtls[1], 4), " ", round(qtls, 4), "]\n", sep = ""))
## find the posterior mean, variance and std. deviation
## using the output from the previous example
post.mean = mean(results)
post.var = var(results)
post.sd = sd(results)
# calculate an approximate 95% credible region using the posterior mean and
# std. deviation
lb = post.mean - qnorm(0.975) * post.sd
ub = post.mean + qnorm(0.975) * post.sd
cat(paste("Approximate 95% credible interval : ["
, round(lb, 4), " ", round(ub, 4), "]\n", sep = ""))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.