par(mfrow=c(1,1)) integrand=function(u,beta){ return(dpois(8,exp(beta+u))*dnorm(u)) } betas=c(-10:40)/10 quadlike=rep(0,length(betas)) for (i in 1:length(betas)){ quadlike[i]=integrate(integrand,-10,10,beta=betas[i])$value } plot(betas,quadlike,type="l") #exercise 1 #be careful with errors in LSMC LSMC=function(y,beta,M){ u=rnorm(1) lambda=exp(beta+u) return(list(mcestimate=mean(dpois(y,lambda*u)),mcstderr=sqrt(var(dpois(y,lambda)/M)))) } SMClike=rep(0,length(betas)) SMCerr=rep(0,length(betas)) M=1000 for (i in 1:length(betas)){ temp=LSMC(8,betas[i],M) SMClike[i]=temp$mcestimate SMCerr[i]=temp$mcstderr } lines(betas,SMClike,lty=2) #now more observations y=c(8,18,5,7,10,9,9,6,7,10) u=seq(-1,5,len=100) densityprod=function(lambda){ #we standardise with poisson density at lambda=mean(y) for numerical stability return(exp( sum(y)*(log(lambda)-log(mean(y))) - length(y)*(lambda - mean(y)))) } integrand=function(u,beta){ return(densityprod(exp(u+beta))*dnorm(u)) } betas=c(-10:40)/10 quadlike=rep(0,length(betas)) low=-10 up=10 for (i in 1:length(betas)){ quadlike[i]=integrate(integrand,low,up,beta=betas[i])$value } plot(betas,quadlike,type="l") #careful with errors in LSMC LSMC=function(y,beta,M){ u=rnorm(M) lambda=exp(beta+u) return(list(mcestimate=mean(densityprod(u)),mcstderr=sqrt(var(densityprod(lambda)/M)))) } SMClike=rep(0,length(betas)) SMCerr=rep(0,length(betas)) M=1000 for (i in 1:length(betas)){ temp=LSMC(8,betas[i],M) SMClike[i]=temp$mcestimate SMCerr[i]=temp$mcstderr } plot(betas,quadlike,type="l",lty=1) lines(betas,SMClike,lty=2) #Exercise 2 y=c(8,18,5,7,10,9,9,6,7,10) minusg=function(u,beta){ #again we standardise with poisson density at lambda=mean(y) for numerical stability return(-sum(y)*((beta+u)-log(mean(y)))+length(y)*(exp(beta+u)-mean(y))-log(dnorm(u))) } impsampquotient=function(v,beta){ return(exp( -minusg(v,beta)-log( dt((v-mu)/sigma,nu)/sigma)) ) } #plot numerator and denominator in importance sampling quotiont #something wrong in the next 6 lines... nu=10 beta=2 mu=nlm(minusg,0,beta=beta)$estimate sigma=1/sqrt(length(y)*exp(beta+mu)+1) #plot imp.sampling quotient and imp. sampling density plot(u,exp(-minusg(u,beta))/(dt((u-mu)/sigma,nu)/sigma),type="l",lty=1) lines(u,dchisq((u-mu)/sigma,nu)/sigma,lty=2) V=rt(1000,nu) impsamplelike=rep(0,length(betas)) for (i in 1:length(betas)){ mu=nlm(minusg,0,beta=betas[i])$estimate sigma=1/sqrt(length(y)*exp(beta+mu)+1) v=mu+sigma*V impsamplelike[i]=mean(impsampquotient(v,beta)) } #plot importance sampling estimate of likelihood with the likelihood computed using numerical quadrature in Exercise 1 e). plot(betas,log(quadlike),type="l",lty=1) lines(betas,log(impsamplelike),lty=2) #result if we use the same importance distribution (corresponding to beta=2) for all other beta. V=rt(1000,nu) impsamplelike=rep(0,length(betas)) for (i in 1:length(betas)){ mu=nlm(minusg,0,beta=4)$estimate sigma=1/sqrt(length(y)*exp(4+mu)+1) v=mu+sigma*V impsamplelike[i]=mean(impsampquotient(v,betas[i])) } #plot importance sampling estimate of likelihood with the likelihood computed using numerical quadrature in Exercise 1 e). plot(betas,log(quadlike),type="l",lty=1) lines(betas,log(impsamplelike),lty=2)