# blood1.r # # BLOOD1(MODEL,MAXITER,S) # Examining the Coagulation of Blood example on p 284 and p385 of Gelman # et # al., which has joint posterior density: # theta,mu,log(sigma),log(tao)|y ~ prod(j=1:4)N(thetaj|mu,tao^2)* # tao*prod(j=1:4)prod(i=1:n_j)N(y_ij|thetaj,sigma^2) # using: # 1. mode finding algorithm (if model==1) # 2. Gibb's Sampler (if model==2) and number of Gibb's seq=1 - If # >1, see bloodk.r # INPUTS: # MODEL - 1=Mode Finding, 2=Gibb's Sampler # MAXITER - # of iterations to run - (I like 20-100) # S - Used for Model 2 - Make inference from the S_th spot in the sequence to the end (I like .5) # blood1<-function(model,maxiter,S) { # Define algorithm parameters #model<-1 #maxiter<-5 # MODEL specific params # 1. MODE FINDING ... # 2. GIBB'S SAMPLING #S<-2/3 # Make inference from the S_th spot in the Gibb's sequence to the end # define the data, y_ij a<-c(62,60,63,59,NA,NA,NA,NA) b<-c(63,67,71,64,65,66,NA,NA) c<-c(68,66,71,67,68,68,NA,NA) d<-c(56,62,60,61,63,64,63,59) m<-8 n<-4 lenvec<-c(length(a[(!is.na(a))]),length(b[(!is.na(b))]),length(c[(!is.na(c))]),length(d[(!is.na(d))])) N<-sum(lenvec) # total samples across n experiments lenmat<-t(matrix(c(4,6,6,8),n,m)) x<-matrix(c(a,b,c,d),m,n) ymj<-apply(x,2,mean,na.rm=T) # ********************************************* # Use mode finding algorithm to simulate the jp # Get first estimates for thetaj, mu, sigma and tao if (model==1) { thetaj<-ymj # initial estimate for theta thetam<-t(matrix(thetaj,n,m)) sigmaj2<-apply((x-thetam)^2/(lenmat-1),2,sum,na.rm=T) sigma2<-mean(sigmaj2) sigma<-sqrt(sigma2) mu<-mean(thetaj) tao2<-var(thetaj) tao<-sqrt(tao2) thetahat<-matrix(0,n,maxiter+1) thetahat[1:n,1]<-thetaj muhat<-seq(0,0,length=maxiter+1) muhat[1]<-mu sigmahat<-seq(0,0,length=maxiter+1) sigmahat[1]<-sigma taohat<-seq(0,0,length=maxiter+1) taohat[1]<-tao jp<-seq(0,0,length=maxiter+1) jp[1]<-tao*prod(dnorm(thetaj,mu,tao))*prod(apply(dnorm(x,thetam,sigma),2,prod,na.rm=T)) for (i in 2:(maxiter+1)) { # Get conditional modes for each thetaj # (see 9.11 and 9.12 on p 286) thetahat[1:n,i]<-(1/tao2*mu + lenvec/sigma2*ymj)/(1/tao2 + lenvec/sigma2) thetaj<-thetahat[1:n,i] thetam<-t(matrix(thetaj,n,m)) # Get conditional modes for mu (9.15) muhat[i]<-1/n*sum(thetaj) mu<-muhat[i] # Get conditional modes for sigma (9.17) sigmahat[i]<-sqrt(1/N*sum(apply((x-thetam)^2,2,sum,na.rm=T))) sigma<-sigmahat[i] sigma2<-sigma^2 # Get conditional modes for tao taohat[i]<-sqrt(1/(n-1)*sum((thetaj-mu)^2)) tao<-taohat[i] tao2<-tao^2 # Get log(p(theta,mu,log(sigma),log(tao)|y)) = # tao*prod(j=1:4)prod(i=1:n_j)N(y_ij|thetaj,sigma^2) jp[i]<-tao*prod(dnorm(thetaj,mu,tao))*prod(apply(dnorm(x,thetam,sigma),2,prod,na.rm=T)) } print("Using mode finding algorithm to simulate p(theta,mu,sigma,tao|y)") } # ********************************************* # Use Gibbs Sampler to simulate the jp # Get first estimates for thetaj, mu, sigma and tao if (model==2) { thetaj<-ymj # initial estimate for theta thetam<-t(matrix(thetaj,n,m)) sigmaj2<-apply((x-thetam)^2/(lenmat-1),2,sum,na.rm=T) sigma2<-mean(sigmaj2) sigma<-sqrt(sigma2) mu<-mean(thetaj) tao2<-var(thetaj) tao<-sqrt(tao2) thetahat<-matrix(0,n,maxiter+1) thetahat[1:n,1]<-thetaj muhat<-seq(0,0,length=maxiter+1) muhat[1]<-mu sigmahat<-seq(0,0,length=maxiter+1) sigmahat[1]<-sigma taohat<-seq(0,0,length=maxiter+1) taohat[1]<-tao jp<-seq(0,0,length=maxiter+1) jp[1]<-tao*prod(dnorm(thetaj,mu,tao))*prod(apply(dnorm(x,thetam,sigma),2,prod,n)) for (i in 2:(maxiter+1)) { # Get theta from N(Mtheta,Vtheta) (9.11) Mtheta<-(1/tao2*mu + lenvec/sigma2*ymj)/(1/tao2 + lenvec/sigma2) Vtheta<-1/(1/tao2 + lenvec/sigma2) thetahat[1:n,i]<-rnorm(n,Mtheta,sqrt(Vtheta)) thetaj<-thetahat[1:n,i] thetam<-t(matrix(thetaj,n,m)) # Get mu from N(Mmu,Vmu) (9.14) Mmu<-1/n*sum(thetaj) Vmu<-tao2/n muhat[i]<-rnorm(1,Mmu,sqrt(Vmu)) mu<-muhat[i] # Get sigma from scaled Inv-ci^2(df=Dfsigma,scale=Ssigma) (9.16) Dfsigma<-N Ssigma2<-1/N*sum(apply((x-thetam)^2,2,sum,na.rm=T)) # scaled Inv-chi^2=Df*scale^2/chisq sigmahat[i]<-sqrt(Dfsigma*Ssigma2/rchisq(1,Dfsigma)) # scaled Inv-chi^2=Inv-gamma #sigmahat[i]<-1/sqrt(rgamma(1,Dfsigma/2,Dfsigma/2*Ssigma2)) sigma<-sigmahat[i] sigma2<-sigma^2 # Get tao from scaled Inv-ci^2(df=Dftao,scale=Stao) (9.18) Dftao<-n-1 Stao2<-1/(n-1)*sum((thetaj-mu)^2) # scaled Inv-chi^2=scale^2/chisq taohat[i]<-sqrt(Dftao*Stao2/rchisq(1,Dftao)) # scaled Inv-chi^2=Inv-gamma #taohat[i]<-1/sqrt(rgamma(1,Dftao/2,Dftao/2*Stao2)) tao<-taohat[i] tao2<-tao^2 # Get log(p(theta,mu,log(sigma),log(tao)|y)) = # tao*prod(j=1:4)prod(i=1:n_j)N(y_ij|thetaj,sigma^2) (p 285) jp[i]<-tao*prod(dnorm(thetaj,mu,tao))*prod(apply(dnorm(x,thetam,sigma),2,prod,na.rm=T)) } print("Using Gibb's Sampler to simulate p(theta,mu,sigma,tao|y)") } # Compute potential scale reduction, sqrt(R) #psimean #B<-N/(n-1)*sum # Show numerical results of the method tot<-matrix(0,n+4,maxiter+1) tot[1:n,1:(maxiter+1)]<-thetahat tot[n+1,1:(maxiter+1)]<-muhat tot[n+2,1:(maxiter+1)]<-sigmahat tot[n+3,1:(maxiter+1)]<-taohat tot[n+4,1:(maxiter+1)]<-log(jp) rownames(tot)<-c("theta1","theta2","theta3","theta4","mu","sigma","tao","log(jp)") #dimnames(tot)<-list(c("theta1","theta2","theta3","theta4","mu","sigma","tao","log(jp)"),c("crude","1st","2nd","3rd","4th")) print(tot) if (model==2) { print("Posterior Inference for individual parameters") probs<-c(.025,.25,.5,.75,.975) suminf<-matrix(0,n+4,6) for (i in 1:n) { suminf[i,1:5]<-quantile(thetahat[i,floor((maxiter+1)*S):(maxiter+1)],probs) } suminf[n+1,1:5]<-quantile(muhat[floor((maxiter+1)*S):(maxiter+1)],probs) suminf[n+2,1:5]<-quantile(sigmahat[floor((maxiter+1)*S):(maxiter+1)],probs) suminf[n+3,1:5]<-quantile(taohat[floor((maxiter+1)*S):(maxiter+1)],probs) dimnames(suminf)<-list(c("theta1","theta2","theta3","theta4","mu","sigma","tao","log(jp)"),c(2.5,25,50,75,97.5,"sqrt(R)")) print(suminf) } }