yx <-read.table("c:/jorge/business/bayes/trans3_1_4g1.txt") y<-as.matrix(yx[,1]) x<-as.matrix(yx[,2:16]) nobs<-nrow(x) nvar<-ncol(x) sig2<-0.05 rstar<-0.85 v<-matrix(data=0.05, nrow=nobs, ncol=1) ones<-matrix(data=1, nrow=nobs, ncol=1) effic<-matrix(data=NA, nrow=nobs, ncol=1) samp<-matrix(data=NA, nrow=10000, ncol=nobs+nvar+2) #(1)STARTING BIG SAMPLING LOOP for(i in 1:11000) { ystar<- y - v betahat<-as.matrix( solve(t(x) %*% x) %*% t(x) %*% ystar) covb<-sig2*solve(t(x)%*%x) beta <- as.matrix(mvrnorm(n=1, betahat, covb, tol=1e-6, empirical=FALSE)) #(2)SAMPLE THE PRECISION (h) FROM AN CHI-SQUARE #CONDITIONAL ON BETA AND ESTIMATE THE VARIANCE e<-ystar-x%*%beta sse<-t(e)%*%e h<- rgamma(1,shape=(nobs-2)/2,scale=sse/2) sig2<-1/h #(3) SAMPLE LAMBDA, THE MEAN OF THE INFFICIENCY #ERROR FROM ITS CONDITIONAL df1<- nobs+1 df2<-1/(t(v)%*%ones-log(rstar)) lambdainv<-rgamma(1,shape=df1,scale=df2) lambda<-1/lambdainv #(4)FIND THE CONDITIONAL MEAN OF THE TRUNCATED #NORMAL, I.E., SAMPLE V FROM ITS CONDITIONAL #SAMPLE FROM THE TRUNCATED NORMAL meanv<-y-x%*%beta-(sig2/lambda)*ones varv<-sig2*ones b<-meanv/sqrt(varv) p<-matrix(data=NA, ncol=1, nrow=nobs) for(j in 1:nobs) { c1<- b[j] t4<-0.450 if(c1>t4) { z<- -log(runif(1)/c1) while(runif(1)>exp(-0.5*z*z)) { z<- -log(runif(1))/c1 } p[j]<-c1+z } else { p[j]<- runif(1) while(p[j]1000) { samp[i-1000, ]<-c( t(effic),t(beta),sig2,lambda) } } mean(samp[,86]) hist(samp[,86]) mean(samp[,10]) hist(samp[,10])