[R] Bayesian inference: Poisson distribution with normal (!) prior

[Carsten Steinhoff]

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)