B.mat  <- matrix(summary(norr.tob)$coef,nrow=1)     # STARTING POINTS
s.sq   <- matrix(rep(10,6), nrow=1)                 # PRIOR on SIGMA^2 OF BETA
y.star <- Y                                         # LATENT DATA START
for (i in mc.start:mc.stop)  {
    for (j in 1:length(Y)) {
        if (Y[j] == 0)  y.star[j] <- rtnorm(1,X[j,]%*%B.mat[i,],s.sq[i,],lower=-Inf,upper=0)
    }
    delta <- t(y.star - X %*% B.mat[i,]) %*% (y.star - X %*% B.mat[i,])
    s.temp <- rep(Inf,ncol(s.sq)); while(!is.finite(sum(s.temp)))
                    s.temp <- 1/rgamma(ncol(s.sq), (gamma0+length(Y))/2, (gamma1+delta)/2 )
    s.sq <- rbind(s.sq,s.temp)
    B.hat <- solve( B0 + t(X)%*%X ) %*% (B0*beta0 + t(X)%*%y.star)
    B.var <-  make.symmetric( solve( s.sq[(i+1)]^(-1)*B0 + s.sq[(i+1)]^(-1)*t(X)%*%X ) )
    B.mat <- rbind( B.mat, rmultinorm(1,B.hat,B.var,tol=1e-06) )
    }
}