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]