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

From: Carsten Steinhoff <carsten.steinhoff_at_gmx.de>
Date: Fri 26 Jan 2007 - 13:56:19 GMT


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)



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 01:04:31 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.