bb=function(nobs=10^2,ncomp=5,npart=10^3, lam=c(10,30,110,150,220,250,310),p=c(.5,.5,.1,.1,.1,.05,.05), nsim=250,obyo=F,wa1=T,seedv=NULL,sobs=c(100,1000,5000,10000)){ library(combinat) # Generating a sample from a mixture of ncomp Poissons # p=c(.5,.5,.1,.1,.1,.05,.05) # lam=c(10,30,110,150,220,250,310) # nobs=10^2 # ncomp=4 # npart=10^3 #---------------------------------------------------------------- # Poisson mixture likelihood like=function(pars,curobs=nobs,obs=obs){ #---------------------------------------------------------------- dens=pars[1]*dpois(obs[1:curobs],pars[ncomp+1]) for (v in 2:ncomp) dens=dens+pars[v]*dpois(obs[1:curobs],pars[ncomp+v]) sum(log(dens)) } #---------------------------------------------------------------- #---------------------------------------------------------------- # Regular MCMC tuned for sequential sampling MCMC=function(start,curobs=nobs,obs=obs){ #---------------------------------------------------------------- zobs=obs[1:curobs] samc=matrix(0,ncol=2*ncomp,nrow=npart) samc[1,]=start dens=matrix(0,nrow=curobs,ncol=ncomp) siz=sumz=matrix(0,nrow=npart,ncol=ncomp) suff=matrix(0,nrow=2*ncomp,ncol=npart) for (t in 2:npart){ #completion dens[,1]=samc[t-1,1]*dpois(zobs,samc[t-1,ncomp+1]) for (v in 2:ncomp) dens[,v]=dens[,v-1]+samc[t-1,v]*dpois(zobs,samc[t-1,ncomp+v]) unive=runif(curobs)*dens[,ncomp] labl=rep(1,curobs) for (v in 1:(ncomp-1)) labl=labl+(unive>dens[,v]) for (v in 1:ncomp){ siz[t-1,v]=sum(labl==v) sumz[t-1,v]=sum(zobs*(labl==v)) } #propagation for (v in 1:ncomp){ samc[t,v]=rgamma(1,gam[v]+siz[t-1,v]) samc[t,ncomp+v]=rgamma(1,alpha[v]+sumz[t-1,v])/(beta[v]+siz[t-1,v]) } samc[t,1:ncomp]=samc[t,1:ncomp]/sum(samc[t,1:ncomp]) #reorderin samc[t,1:ncomp]=samc[t,order(samc[t,ncomp+(1:ncomp)])] samc[t,ncomp+(1:ncomp)]=sort(samc[t,ncomp+(1:ncomp)]) } #completion dens[,1]=samc[t-1,1]*dpois(zobs,samc[t-1,ncomp+1]) for (v in 2:ncomp) dens[,v]=dens[,v-1]+samc[t-1,v]*dpois(zobs,samc[t-1,ncomp+v]) unive=runif(curobs)*dens[,ncomp] labl=rep(1,curobs) for (v in 1:(ncomp-1)) labl=labl+(unive>dens[,v]) for (v in 1:ncomp){ siz[npart,v]=sum(labl==v) sumz[npart,v]=sum(zobs*(labl==v)) } for (v in 1:ncomp){ suff[v,]=siz[,v];suff[ncomp+v,]=sumz[,v]} list(paras=samc,hypsiz=siz,hypsum=sumz,suff=suff) } #---------------------------------------------------------------- #---------------------------------------------------------------- # Chib's approximation Chib=function(parm,hyper,curobs=nobs,obs=obs){ #---------------------------------------------------------------- logliu=apply(parm,1,like,curobs=curobs,obs=obs) ovid=max(logliu) logliu=order(logliu)[npart] pstar=parm[logliu,1:ncomp] lastar=parm[logliu,ncomp+(1:ncomp)] aumone=matrix(0,ncol=factorial(ncomp),nrow=npart) for (n in 1:factorial(ncomp)){ origo=permn(1:ncomp)[[n]] denome=rep(lgamma(sum(gam)+curobs),npart) for (v in (1:ncomp)) denome=denome+ (gam[v]+hyper[v,]-1)*log(pstar[origo[v]])- lgamma(gam[v]+hyper[v,])+(alpha[v]+hyper[ncomp+v,]-1)* log(lastar[origo[v]])-(beta[v]+hyper[v,])*lastar[origo[v]]- lgamma(alpha[v]+hyper[ncomp+v,])+(alpha[v]+hyper[ncomp+v,]-1)* log(beta[v]+hyper[v,]) aumone[,n]=denome } ovid+lgamma(sum(gam))+sum((gam-1)*log(pstar))-sum(lgamma(gam))+ sum(log(beta)+dgamma(lastar*beta,alpha,log=TRUE))- log(mean(exp(as.vector(aumone)))) } #---------------------------------------------------------------- #********** # MAIN #********** if(length(seedv)weight[newindices,v]) for (v in 1:ncomp){ suff[v,newz==v]=suff[v,newz==v]+1 suff[ncomp+v,newz==v]=suff[ncomp+v,newz==v]+obs[t] } # generating parameters psim=lsim=matrix(0,ncol=ncomp,nrow=npart) for (v in 1:ncomp){ psim[,v]=rgamma(npart,gam[v]+suff[v,]) lsim[,v]=rgamma(npart,alpha[v]+suff[ncomp+v,])/(beta[v]+suff[v,]) } psim=psim/apply(psim,1,sum) if(obyo||sum(t==sobs)>0){ # Chib's approximation based on PL evid[t,2]=Chib(parm=cbind(psim,lsim),hype=suff,curobs=t,obs) # MCMC comparison, same number of iterations mcmres=MCMC(stard,curobs=t,obs) stard=mcmres\$paras[npart,] # Chib's approximation based on MCMC evid[t,3]=Chib(parm=mcmres\$par,hype=mcmres\$suff,curobs=t,obs) } # end if (obyo...) } # end t nfall=paste("all.",seedv[isim],".dat",sep="") fall=file(nfall,"w") if(wa1){ # write all values of evid[,1] write(t(evid),ncol=3,file=fall) }else if(!is.null(sobs)){ #write only values at t=sobs write(t(evid[sobs,]),ncol=3,file=fall) } close(fall) if(wa1) system(paste("gzip",nfall,sep=" ")) } # end isim } # end bb(){