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]<c1)
			{
				p[j]<- runif(1)
			}
		}
	}
	v<-meanv+p*sqrt(varv)
	effic<- exp(-v)

	#DROP INITIAL 1000 POINTS AND PUT REST IN SAMP
	if(i>1000)
	{
		samp[i-1000, ]<-c( t(effic),t(beta),sig2,lambda)
	}
}
	
mean(samp[,86])
hist(samp[,86])
mean(samp[,10])
hist(samp[,10])