poisgcp {Bolstad} | R Documentation |
Poisson sampling with a general continuous prior
Description
Evaluates and plots the posterior density for mu, the mean rate
of occurance of an event or objects, with Poisson sampling and a
general continuous prior on mu
Usage
poisgcp(y, density = "normal", params=c(0,1), n.mu = 100
,mu = NULL, mu.prior = NULL,
print.sum.stat=FALSE, alpha=0.05, ret=FALSE)
Arguments
y |
A random sample of one or more observations from a
Poisson distribution |
density |
may be one of "gamma", "normal", or "user" |
params |
if density is one of the parameteric forms then
then a vector of parameters must be supplied.
gamma: a0,b0
normal: mean,sd |
n.mu |
the number of possible mu values in the prior. This
number must be greater than or equal to 100. It is ignored
when density="user". |
mu |
a vector of possibilities for the mean of a Poisson distribution. This must be set if density="user". |
mu.prior |
the associated prior density. This must be set if density="user". |
print.sum.stat |
if set to TRUE then the posterior mean,
posterior variance, and a credible interval for the mean are
printed. The width of the credible interval is controlled by
the parameter alpha. |
alpha |
The width of the credible interval is controlled by
the parameter alpha. |
ret |
this arguement is deprecated. |
Value
A list will be returned with the following components:
mu |
the vector of possible mu values used in the prior |
mu.prior |
the associated probability mass for the values in mu |
likelihood |
the scaled likelihood function for
mu given y |
posterior |
the posterior probability of mu given y |
See Also
poisdp
poisgamp
Examples
## Our data is random sample is 3, 4, 3, 0, 1. We will try a normal
## prior with a mean of 2 and a standard deviation of 0.5.
y = c(3,4,3,0,1)
poisgcp(y,density="normal",params=c(2,0.5))
## The same data as above, but with a gamma(6,8) prior
y = c(3,4,3,0,1)
poisgcp(y,density="gamma",params=c(6,8))
## The same data as above, but a user specified continuous prior.
## We will use print.sum.stat to get a 99% credible interval for mu.
y = c(3,4,3,0,1)
mu = seq(0,8,by=0.001)
mu.prior = c(seq(0,2,by=0.001),rep(2,1999),seq(2,0,by=-0.0005))/10
poisgcp(y,"user",mu=mu,mu.prior=mu.prior,print.sum.stat=TRUE,alpha=0.01)
## find the posterior CDF using the results from the previous example
## and Simpson's rule
results = poisgcp(y,"user",mu=mu,mu.prior=mu.prior)
cdf = sintegral(mu,results$posterior,n.pts=length(mu),ret=TRUE)
plot(cdf,type="l",xlab=expression(mu[0])
,ylab=expression(Pr(mu<=mu[0])))
## use the cdf to find the 95% credible region.
lcb = cdf$x[with(cdf,which.max(x[y<=0.025]))]
ucb = cdf$x[with(cdf,which.max(x[y<=0.975]))]
cat(paste("Approximate 95% credible interval : ["
,round(lcb,4)," ",round(ucb,4),"]\n",sep=""))
## find the posterior mean, variance and std. deviation
## using Simpson's rule and the output from the previous example
dens = mu*results$posterior # calculate mu*f(mu | x, n)
post.mean = sintegral(mu,dens)
dens = (mu-post.mean)^2*results$posterior
post.var = sintegral(mu,dens)
post.sd = sqrt(post.var)
# 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=""))
[Package
Bolstad version 0.2-17
Index]