poisgcp {Bolstad}

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 = c("normal", "gamma", "user"),
  params = c(0, 1),
  n.mu = 100,
  mu = NULL,
  mu.prior = NULL,
  print.sum.stat = FALSE,
  alpha = 0.05,
  ...
)

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

either a vector of possibilities for the mean of a Poisson distribution, or a range (a vector of length 2) of values. This must be set if density = "user". If mu is a range, then n.mu will be used to decide how many points to discretise this range over.

mu.prior

either a vector containing y values correspoding to the values in mu, or a function. This is used to specifiy the prior f(\mu). So mu.prior can be a vector containing f(\mu_i) for every \mu_i, or a funtion. 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.

...

additional arguments that are passed to Bolstad.control

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. Note that the syntax of sintegral has changed.
results = poisgcp(y,"user",mu=mu,mu.prior=mu.prior)
cdf = sintegral(mu,results$posterior,n.pts=length(mu))$cdf
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)$value

dens = (mu-post.mean)^2*results$posterior
post.var = sintegral(mu,dens)$value
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=""))

# NOTE: All the examples given above can now be done trivially in this package

## find the posterior CDF using the results from the previous example
results = poisgcp(y,"user",mu=mu,mu.prior=mu.prior)
cdf = cdf(results)
curve(cdf,type="l",xlab=expression(mu[0])
	,ylab=expression(Pr(mu<=mu[0])))

## use the quantile function to find the 95% credible region.
ci = quantile(results, c(0.025, 0.975))
cat(paste0("Approximate 95% credible interval : ["
	,round(ci[1],4)," ",round(ci[2],4),"]\n"))

## 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
ci = post.mean + c(-1, 1) * qnorm(0.975) * post.sd

cat(paste("Approximate 95% credible interval : ["
	,round(ci[1],4)," ",round(ci[2],4),"]\n",sep=""))

## Example 10.1 Dianna's prior
# Firstly we need to write a function that replicates Diana's prior
f = function(mu){
   result = rep(0, length(mu))
   result[mu >=0 & mu <=2] = mu[mu >=0 & mu <=2]
   result[mu >=2 & mu <=4] = 2
   result[mu >=4 & mu <=8] = 4 - 0.5 * mu[mu >=4 & mu <=8]
   
   ## we don't need to scale so the prior integrates to one, 
   ## but it makes the results nicer to see
   
   A = 2 + 4 + 4
   result = result / A
   
   return(result)
 }
 
 results = poisgcp(y, mu = c(0, 10), mu.prior = f)