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

[Vincent Goulet]

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