## 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

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)
,ylab=expression(Pr(mu<=mu)))

## 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)
,ylab=expression(Pr(mu<=mu)))

## 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,4)," ",round(ci,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,4)," ",round(ci,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)