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

From: Vincent Goulet <vincent.goulet_at_act.ulaval.ca>
Date: Fri 26 Jan 2007 - 15:14:14 GMT

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

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Sat Jan 27 02:22:55 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 26 Jan 2007 - 15:30:29 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.