[R] Bayesian inference: Poisson distribution with normal (!) prior
Vincent Goulet
vincent.goulet at act.ulaval.ca
Fri Jan 26 16:14:14 CET 2007
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
More information about the R-help
mailing list