llogitpostall1 <- function(data,pia,pib,etaa,etab,lama,lamb,lpi,leta,llambda)
   {
    # updated August 29, 2006
     samp     <- nrow(data)
     posprob <- numeric(samp)
     negprob <- numeric(samp)
     lres    <- numeric(samp)
     r       <- data[,1]
     neg     <- data[,2]
     pos     <- data[,3]
     pi      <-  exp(lpi)/(1 + exp(lpi) )
     eta     <-  exp(leta)/(1 + exp(leta) )
     lambda  <-  exp(llambda)/(1 + exp(llambda) )

     negprob <-  ( exp(llambda)/(1+exp(llambda))*(1+exp(lpi))^(-r)+(1+exp(leta))^(-1)*(1 -(1 +exp(lpi))^(-r)))
     posprob <-  (exp(leta)/(1+exp(leta)))*(1 -(1 +exp(lpi))^(-r))+(1+exp(llambda))^(-1)*(1+exp(lpi))^(-r)
     jacobian <- (exp(llambda)/(1+exp(llambda))^2)*(exp(leta)/(1+exp(leta))^2)*(exp(lpi)/(1+exp(lpi))^2)
     lres <- neg*log(negprob) + pos*log(posprob) 

     prior   <- (pia-1)*log(pi) + (pib-1)*log((1-pi)) +
                (etaa-1)*log(eta) + (etab-1)*log((1-eta)) +
                (lama-1)*log(lambda) + (lamb-1)*log((1-lambda))

     ltot    <- sum(lres) +prior +log(jacobian)
    # print("r, neg, pos ") ; print(cbind(r,neg,pos))
    # print("posprob, negprob ") ; print(cbind(posprob,negprob))
    # print("lres ") ; print(lres)
    # print("prior ") ; print(prior)
    # print("ltot ") ; print(ltot)
     ltot
   }

convplt1 <- function(y,wl)
   {
     num <- length(y)/wl
     out25 <- numeric(num)
     out50 <- numeric(num)
     out75 <- numeric(num)
     dim(y) <- c(wl,num)
     out25 <- apply(y,2,quantile,probs=.25)
     out50 <- apply(y,2,quantile,probs=.5)
     out75 <- apply(y,2,quantile,probs=.75)
     plot(c(1:num),out25,type="b",pch="1",xlab="iterations",ylab="quantiles",
           ylim=c(min(out25),max(out75)))
     lines(c(1:num),out50,type="b",pch="2")
     lines(c(1:num),out75,type="b",pch="3")
   }


poolprevsample1 <- function(data,pia,pib,etaa,etab,lama,lamb,burnin,finalsamp,
                    std=c(1.,1.,1.),start=c(0.,0.,0.),output=1)
   {
     begin <- proc.time()
     iter  <- burnin +finalsamp
     logitpisamp <- numeric(iter)
     logitetasamp <- numeric(iter)
     logitlambsamp <- numeric(iter)
     accept  <- numeric(iter) 
     logitpisamp[1]  <-  start[1] ; logitetasamp[1]  <- start[2] ;  logitlambsamp[1]  <- start[3]
     for (i in 2:iter)
     {
       lproppi  <- rnorm(1,mean=logitpisamp[i-1],std[1])
       lpropeta  <- rnorm(1,mean=logitetasamp[i-1],std[2])
       lproplamb  <- rnorm(1,mean=logitlambsamp[i-1],std[3])
       rdens <- exp(llogitpostall1(data,pia,pib,etaa,etab,lama,lamb,lproppi,lpropeta,lproplamb) -
            llogitpostall1(data,pia,pib,etaa,etab,lama,lamb,logitpisamp[i-1],logitetasamp[i-1],logitlambsamp[i-1]) )
       u  <- runif(1)
       if ( u < min(rdens,1) )
          {
            logitpisamp[i]  <-  lproppi ; logitetasamp[i]  <- lpropeta ;  
            logitlambsamp[i]  <- lproplamb ; accept[i] <- 1
          }     
       if ( u >= min(rdens,1) )
          {
            logitpisamp[i]  <- logitpisamp[i-1] ; logitetasamp[i]  <- logitetasamp[i-1] 
            logitlambsamp[i]  <- logitlambsamp[i-1] ; accept[i] <- 0
          }   
     }
     cutoff <- burnin
     lpifinal  <- logitpisamp[cutoff:iter]
     letafinal  <- logitetasamp[cutoff:iter]
     llambfinal  <- logitlambsamp[cutoff:iter]
     pifinal <- exp(lpifinal)/(1 + exp(lpifinal) )
     etafinal <- exp(letafinal)/(1 + exp(letafinal) )
     lambfinal <- exp(llambfinal)/(1 + exp(llambfinal) )
     # print(c(length(etasamp),length(lambsamp),length(etafinal),length(lambfinal)))
     qpi <- quantile(pifinal,c(.01,.025,.05,.1,.25,.5,.75,.9,.95,.975,.99))
     print("Quantiles for pi (Prevalence)") ; print(qpi)
     qeta <- quantile(etafinal,c(.01,.025,.05,.1,.25,.5,.75,.9,.95,.975,.99))
     print("Quantiles for eta (Sensitivity)") ; print(qeta)
     qlamb <- quantile(lambfinal,c(.01,.025,.05,.1,.25,.5,.75,.9,.95,.975,.99))
     print("Quantiles for lambda (Specificity)") ; print(qlamb)
     if (output == 2)
       { 
          ckcorrpi <- acf(pifinal)
          windows()
          print(c("For pifinal, autocorrelations = ",ckcorrpi$acf[2:11,,]),digits=5)
          convplt1(logitpisamp,1000)
          windows()
          pairs(cbind(lpifinal,letafinal,llambfinal))
          windows()
          pairs(cbind(pifinal,etafinal,lambfinal))
          windows()
          hist(etafinal)
          windows()
          hist(lambfinal)
       }
     hist(pifinal) 
     print("Mean and SD for pi (Prevalence)") ; print(c(mean(pifinal),sd(pifinal)))
     print("Mean and SD for eta (Sensitivity)") ; print(c(mean(etafinal),sd(etafinal)))
     print("Mean and SD for lambda (Specificity)") ; print(c(mean(lambfinal),sd(lambfinal)))
     end <- proc.time()
     print("Proportion of accepted proposals is "); print( sum(accept)/iter) 
     print("total time in minutes is "); print( (end[3]-begin[3])/60) 
     ans <- list(prevquantiles=qpi,prevmean=mean(pifinal),prevsd=sd(pifinal),
                 sensquantiles=qeta,sensmean=mean(etafinal),senssd=sd(etafinal),
                 specquantiles=qlamb,specmean=mean(lambfinal),specsd=sd(lambfinal))
   }