gibbs_ function(X, Y, rows, cols, num.iter = 1000, 
    opt = 0, ptab.len = rows,	
    pta = rep(1.1, length(rows)), ptb = rep(.1,length(rows)), 
    prior.var  = 1000000, g.diag = rep(1.0,nrow(X)),
    burnin = 100)
{
    if (opt != 0 && opt != 1 ) stop("\nopt can only be 0, 1")
    num.err_ length(ptab.len)
    if (num.err != length(pta) || num.err != length(ptb) )
        stop("\npta, ptb, and ptab.len have different lengths")
    rX_ dim(X)[1]; cX_ dim(X)[2]
    pri.var.len_ length(prior.var)
    c.err.len_ sum(ptab.len)
    if (pri.var.len != 1 && (pri.var.len + c.err.len) != rX) 
        stop("\ntoo many or too few error terms")
    ran.err_ rgamma(num.err,pta)/(ptb+10^(-8))
    err_ rep(1/prior.var, rX-c.err.len)
    for (i in num.err:1) err _ c(rep(ran.err[i],ptab.len[i]),err)
    gamm_ g.diag * err
    ahit_ ptab.len/2+pta
    hit_ lhit_ rep(NA,num.err)
    i _ 1
    while (i <= burnin) {	
        sv   _ svd(t(X * gamm) %*% X)
        Vhit _ sv$u %*% (t(sv$v) * (1/sv$d))
        eThit_ Vhit %*% (t(X) %*% (gamm * Y))
        Thit _ sv$u %*% ((t(sv$v) * sv$d^(-.5))
            %*% as.matrix(rnorm(cX))) + eThit
        resd _ g.diag * (Y - X %*% Thit)^2
        ix1_ 0
        for (j in 1:num.err) {
            ix2_ ptab.len[j] + ix1
            lhit[j]_ ptb[j] + sum(resd[(ix1+1):ix2])/2
        ix2_ ix1
        }
    hit_ rgamma(num.err,ahit)/lhit
    err_ rep(1/prior.var, rX-c.err.len)
    for (j in num.err:1) err_ c(rep(hit[j],ptab.len[j]),err)
    gamm_ g.diag * err
    i _ i+1
    }
    lhit_ matrix(NA,num.iter,num.err)
    eThit_ vThit_ matrix(NA,num.iter,cX)
    if (opt==1) {
        Thit _ matrix(NA,num.iter,cX)
        hit_ matrix(NA,num.iter,num.err)
    i_ 1
    while (i <= num.iter) {
        sv_ svd(t(X * gamm) %*% X)
        Vhit_ sv$u %*% (t(sv$v) * (1/sv$d))
        vThit[i,]_ diag(Vhit)
        eThit[i,]_ Vhit %*% (t(X) %*% (gamm * Y))
        Thit[i,]_ sv$u %*% ((sv$d^(-.5) * t(sv$v))
                    %*% as.matrix(rnorm(cX))) + eThit[i,]
        resd_ g.diag * (Y - X %*% Thit[i,])^2
        ix1_ 0
        for (j in 1:num.err) {
            ix2_ ptab.len[j] + ix1
            lhit[i,j]_ ptb[j] + sum(resd[(ix1+1):ix2])/2
            ix1_ ix2
            }
        hit[i,] _ rgamma(num.err,ahit)/lhit[i,]
        err_ rep(1/prior.var, rX-c.err.len)
        for (j in num.err:1) err_ c(rep(hit[i,j],ptab.len[j]),err)
        gamm_ g.diag * err
        i _ i+1
        }
    eth _ apply(eThit,2,mean)
    vth _ apply(vThit,2,mean) + apply(eThit,2,var)
    ehe _ apply(ahit/t(lhit),1,mean)
    vhe _ apply(ahit/t(lhit)^2,1,mean) + apply(ahit/t(lhit),1,var)
    esig _ apply(t(lhit)/(ahit-1),1,mean)
    vsig _ apply(t(lhit)^2/((ahit-1)^2*(ahit-2)),1,mean) +
    apply(t(lhit)/(ahit-1),1,var) 
    z_ list(
        Y = Y, X = X, num.iter = num.iter, 
        rows = rows, cols = cols, ptab.len = ptab.len,
        mean.theta = eth, var.theta = vth, 
        ehe = ehe, vhe = vhe, 
        mean.sig2 = esig, var.sig2= vsig,
        Ghe = gamm^.5, g.diag = g.diag,
        Gibbs.iter=list(thit=Thit,hit=hit)
        )
    z
    }
    else {
    i _ 1
    while (i <= num.iter) {
        sv_ svd(t(X * gamm) %*% X)
        Vhit_ sv$u %*% (t(sv$v) * (1/sv$d))
        vThit[i,]_ diag(Vhit)
        eThit[i,]_ Vhit %*% (t(X) %*% (gamm * Y))
        Thit_ sv$u %*% ((t(sv$v) * (sv$d^(-.5)))
            %*% matrix(rnorm(cX),byrow=T)) + eThit[i,]
        resd_ Y - X %*% as.matrix(Thit)
        resd_ g.diag * resd^2
        ix1_ 0
        for (j in 1:num.err) {
            ix2_ ptab.len[j] + ix1
            lhit[i,j]_ ptb[j] + sum(resd[(ix1+1):ix2])/2
            ix1_ ix2
        }
        hit_ rgamma(num.err,ahit)/lhit[i,]
        err_ rep(1/prior.var, rX-c.err.len) 
        for (j in num.err:1) err_ c(rep(hit[j],ptab.len[j]),err)
        gamm_ g.diag * err
        i _ i+1
    }
    eThit_ eThit[1:num.iter,]
    vThit_ vThit[1:num.iter,]
    lhit _ lhit[1:num.iter,]
    eth _ apply(eThit,2,mean)
    vth _ apply(vThit,2,mean) + apply(eThit,2,var)
    ehe _ apply(ahit/t(lhit),1,mean)
    vhe _ apply(ahit/t(lhit)^2,1,mean) + apply(ahit/t(lhit),1,var)
    esig _ apply(t(lhit)/(ahit-1),1,mean)
    vsig _ apply(t(lhit)^2/((ahit-1)^2*(ahit-2)),1,mean) +
    apply(t(lhit)/(ahit-1),1,var) 
    z_ list(Y = Y, X = X, num.iter = num.iter, 
        rows = rows, cols = cols, ptab.len = ptab.len,
        mean.theta = eth, var.theta = vth, 
        ehe = ehe, vhe = vhe, 
        mean.sig2 = esig, var.sig2 = vsig,
        Ghe = gamm^.5, g.diag = g.diag,)
    z
    }
}