R help - Jan 2007 - Bayesian inference: Poisson distribution with normal (!) prior

Hello,

for a frequency modelling problem I want to combine expert knowledge with
incoming real-life data (which is not available up to now). The frequency
has to be modelled with a poisson distribution. The parameter lambda has to
be normal distributed (for certain reasons we did not NOT choose gamma
althoug it would make everything easier).

I've started with the subsequent two functions to obtain random numbers for
Lambda after the first observed period. My question is now, how to get the
randoms for the n following periods?

Thanks a lot for your hints! Maybe there is an easier way to do the
necessary calculations...?

Carsten

# Function 1: Posterior for the first observation
test.posterior=function(x,observation,p1,p2)
{
f1=function(x,observation,p1,p2)
dpois(observation,qnorm(pnorm(x,p1,p2),p1,p2))*dnorm(x,p1,p2)
integral=integrate(f1,0,Inf,p1=p1,p2=p2,observation=observation)$value
ausgabe=f1(x,observation,p1=p1,p2=p2)/integral
return(ausgabe)
}

# Function 2: Random numbers for Lambda in the second period
test.posterior.random=function(n,x,length,observation,p1,p2)
{
# n = random numbers to calculate
# x = maximum value for integral calculation
ret=c()
x=seq(0.001,x,length.out=length)
for (i in x)
{
ret=c(ret,integrate(test.posterior,observation=observation,p1=p1,p2=p2,lower
=1,i)$value)
}
ret=abs(ret)
pr=cbind(ret,x)
pr=which(pr[,1]==unique(pr[,1]))
k=approxfun(ret[pr],x[pr])
return(k(runif(n)))
}

# Generate 1000 random numbers
test.posterior.random(1000,.5,1000,1,.2,.05)

Le Vendredi 26 Janvier 2007 08:56, Carsten Steinhoff a
?crit?:
> Hello,
>
> for a frequency modelling problem I want to combine expert knowledge with
> incoming real-life data (which is not available up to now). The frequency
> has to be modelled with a poisson distribution. The parameter lambda has to
> be normal distributed (for certain reasons we did not NOT choose gamma
> althoug it would make everything easier).
>
> I've started with the subsequent two functions to obtain random numbers
for
> Lambda after the first observed period. My question is now, how to get the
> randoms for the n following periods?
>
> Thanks a lot for your hints! Maybe there is an easier way to do the
> necessary calculations...?
>
> Carsten
>
> # Function 1: Posterior for the first observation
> test.posterior=function(x,observation,p1,p2)
> {
> f1=function(x,observation,p1,p2)
> dpois(observation,qnorm(pnorm(x,p1,p2),p1,p2))*dnorm(x,p1,p2)
> integral=integrate(f1,0,Inf,p1=p1,p2=p2,observation=observation)$value
> ausgabe=f1(x,observation,p1=p1,p2=p2)/integral
> return(ausgabe)
> }
>
> # Function 2: Random numbers for Lambda in the second period
> test.posterior.random=function(n,x,length,observation,p1,p2)
> {
> # n = random numbers to calculate
> # x = maximum value for integral calculation
> ret=c()
> x=seq(0.001,x,length.out=length)
> for (i in x)
> {
> ret=c(ret,integrate(test.posterior,observation=observation,p1=p1,p2=p2,lowe
>r =1,i)$value)
> }
> ret=abs(ret)
> pr=cbind(ret,x)
> pr=which(pr[,1]==unique(pr[,1]))
> k=approxfun(ret[pr],x[pr])
> return(k(runif(n)))
> }
>
> # Generate 1000 random numbers
> test.posterior.random(1000,.5,1000,1,.2,.05)
Ok, I'll not question the idea to use a distribution defined on real numbers
as a prior for a strictly positive parameter... Now, if what you want are 
random numbers from the Poisson/normal mixture, this will do:
> rpois(1000, lambda = rnorm(1000, mean, sd))
The reason this works is because the r* functions are vectorized, so here 
rpois() will generate the first number with lambda equal to the first value 
returned by rnorm(), the second number with lambda equal to the second value 
returned by rnorm(), etc.

Am I missing what you want to do?

-- 
  Vincent Goulet, Associate Professor
  ?cole d'actuariat
  Universit? Laval, Qu?bec 
  Vincent.Goulet at act.ulaval.ca   http://vgoulet.act.ulaval.ca