## Constructing a prior for a Poisson mean

I am interested in learning about my son’s cell phone use. Suppose the mean number of text messages that he sends and receives per day is equal to . I will observe the number of text messages for days — we’ll assume that conditional on , are independent Poisson().

How do I construct a prior density on ? Suppose I model my beliefs by use of a gamma density with shape parameter and rate parameter . I want to figure out the values of the prior parameters and .

There are two distinct ways of assessing your prior. One can think about plausible values of the mean number of daily text messages . Alternatively, it may be easier to think about the actual number of text messages — if we assume a gamma prior, then one can show that the predictive density of is given by

Personally, I think it is easier to think about the number of text messages . My best prediction at is 10. After some additional thought, I decide that my standard deviation for is equal to 3.5. One can show that the mean and standard deviation of the predictive density are given by

If one matches my guesses with these expressions, one can show .

To see if these values make any sense, I plot my predictive density. This density is a special case of a negative binomial density where (using the R notation)

size = a, prob = b /(b + 1).

I graph using the following R commands. One can compute P(5 <= y <= 15) = 0.89 which means that on 89% of the days, I believe Steven will send between 5 and 15 messages.

a = 44.4; b = 4.4

plot(0:30,dnbinom(0:30,size=a,prob=b/(b+1)),type=”h”,

xlab=”y”, ylab=”Predictive Prob”, col=”red”, lwd=2)

## Modeling with Cauchy errors

One attractive feature of Bayesian thinking is that one can consider alternative specifications for the sampling density and prior and these alternatives makes one think more carefully about the typical modeling assumptions.

Generally we observe a sample distributed from a sampling density depending on a parameter . If we assign a prior density , then the posterior density is given (up to a proportionality constant) by

Consider the following Cauchy/Cauchy model as an alternative to the usual Normal/Normal model. We let be distributed from a Cauchy density with location and known scale parameter , and we assign a Cauchy prior with location and scale .

Of course, this is not going to be a conjugate analysis — the posterior density has a complicated form. But it is easy to use R to plot and summarize the posterior density.

First we define the posterior density for . The arguments to the function are theta (which could be a vector), the vector of data values, the known scale parameter , and the prior parameters andn . The output is a vector of values of the (unnormalized) posterior density.

posterior=function(theta,data,scale,mu0,tau)

{

f=function(theta) prod(dcauchy(data,theta,scale))

likelihood=sapply(theta,f)

prior=dcauchy(theta, mu0, tau)

return(prior*likelihood)

}

## How many more home runs will Ryan Howard hit?

Here is a basic proportion problem. We are approaching the end of the 2009 baseball season and Ryan Howard of the Philadelphia Phillies currently has 38 home runs in 540 at-bats (opportunities). I project that he will have 85 more at-bats this season. How many additional home runs will he hit?

Let p be the probability that Howard hits a home run in a single at-bat during the 2009 season. We’ll address this prediction problem in three steps. First, we’ll use past data from previous seasons to form a prior distribution for p. Second, we’ll update our beliefs about p using the 2009 data. Last, we’ll predict the number of future at-bats from the posterior predictive distribution. In this process, we’ll illustrate the use of some functions from the LearnBayes package in R.

1. My prior.

In the four previous seasons (2005 through 2008) his home runs/at-bats have been 22/312. 58/581. 47/529, and 48/610. The home run rates are 0.071, 0.100, 0.089, 0.079. Remember that p is the 2009 home run hitting probability. Based on these data, I believe that

(1) P(p < 0.075) = 0.5

(2) P(p < 0.115) = 0.9

In other words, my best guess at p is 0.075 and I’m pretty confident that p is smaller than 0.115. I find the beta(a, b) that matches this information by use of the beta.select function in LearnBayes. I first define the two quantiles as lists — these quantiles are the arguments in beta.select.

> quantile1=list(p=.5,x=.075)

> quantile2=list(p=.9,x=.115)

> beta.select(quantile1, quantile2)

[1] 7.26 85.77

> n=85

> y=0:85

> pred.probs=pbetap(ab,n,y)

## Do graduate students prefer McDonalds?

Today I did an illustration of discrete Bayes for a proportion. I’m interested in the proportion p of graduate students who answer “McDonalds” when asked the question “McDonalds, Wendys, or Burger King?”

I believe p can be one of the five values 0.1, 0.2, 0.3, 0.4, 0.5 and I assign the respective prior weights 1, 2, 5, 10, 5. I define this prior in R:

> p=c(.1,.2,.3,.4,.5)

> prior=c(1,2,5,10,5)

> prior=prior/sum(prior)

> names(prior)=p

> barplot(prior,ylim=c(0,.6),xlab=”p”,main=”PRIOR”)

> plot(s,xlab=”p”,main=”POSTERIOR”)

## Discrete Bayes, part II

In the previous post, I used the R function dexp in my discrete Bayes function. Sometimes you need to make slight changes to the functions in the R stats package for a given problem.

Here is a finite population sampling problem. A small community has 100 voters and you are interested in learning about the number, say M, who will vote for a school levy. You take a random sample of size 20 (without replacement) and 12 are in support of the levy. What have you learned about M?

Here the sampling density is hypergeometric — the probability that will be in favor out of 20 given that M in the population favor the levy — this is given by the hypergeometric probability

There is a function dhyper that computes this probability, but uses a different parametrization. So I write a new function dhyper2 that uses my parametrization.

dhyper2=function(y,M,sample.size=n,pop.size=N)

dhyper(y,M,N-M,n)

## Discrete Bayes

One good way of introducing Bayesian inference is by the use of discrete priors. I recently wrote a simple generic R function that does discrete Bayes for arbitrary choice of a prior and sampling density. I’ll illustrate this function here and in the next posting.

Suppose I’m interested in learning about the rate parameter of an exponential density. We observe a sample from the density defined by

where we are interested in learning about the rate parameter .

In R, I read in the function “discrete.bayes” and a few associated methods by typing

> source(“http://bayes.bgsu.edu/m6480/R/discrete.bayes.functions.R”)