post1mat <- function(data,pia,pib,etaa,etab,lama,lamb,pi,eta,lambda)
   {
     samples <- dim(data)[1]
     posprob <- numeric(samples)
     negprob <- numeric(samples)
     pos <- data[,2]
     neg <- data[,1]
     r   <- data[,3]
     llpart <- numeric(samples)
     for (i in 1:samples)
     {
       posprob[i] <- eta*(1-(1-pi)^r[i]) + (1-lambda)*(1-pi)^r[i]
       negprob[i] <- (1-eta)*(1-(1-pi)^r[i]) + lambda*(1-pi)^r[i]
       llpart[i] <- neg[i]*log(negprob[i]) + pos[i]*log(posprob[i]) 
     }
     lres <- sum(llpart) + (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))
   #  print("Density value =  ") ; print(exp(sum(llpart))) ;
     res  <- exp(lres)
     res
   }


compenummat1a <- function(data,pia,pib,etaa,etab,lama,lamb)
    {
       begin <- proc.time()
       postv <- array(0, dim=c(199,99,99))
       for (pi in seq(.005,.995,by=.005) )
         {  
          for (eta in seq(.01,.99,by=.01) )
              {
              for (lambda in seq(.01,.99,by=.01) )
                   {
                    i <- round((pi -.005)/.005 + 1)
                    j <- round((eta -.01)/.01 + 1)
                    k <- round((lambda -.01)/.01 + 1)
                    postv[i,j,k] <- post1mat(data,pia,pib,etaa,etab,lama,lamb,pi,eta,lambda)
                   } 
               }
           }
       pimarg <- apply(postv, c(1), sum)
       pimargnorm <- pimarg/sum(pimarg)
       print(dim(postv)) ; print(length(pimarg)) 
       print(cumsum(pimargnorm))
       plot(seq(.005,.995,by=.005),pimarg,type="l")
       windows() ; plot(seq(.005,.995,by=.005),pimargnorm,type="l",xlab="Prevalence",ylab="Posterior Density")
      # windows() ; plot(1:80/200,pimargnorm[1:80],type="l")
       quant <- c(.005,.01,.025,.05,.10,.25,.5,.75,.90,.95,.975,.99,.995)
       qres <- approx(x=cumsum(pimargnorm),y=seq(.005,.995,by=.005),xout=quant)
       pimean <- t(seq(.005,.995,by=.005))%*%pimargnorm
       pivar <- t((seq(.005,.995,by=.005)-pimean)^2)%*%pimargnorm
       pistd <- sqrt(pivar)
       print(qres)
       print("Mean and SD for pi") ; print(c(pimean,pistd))
       end <- proc.time()
       print("total time in minutes is "); print( (end[3]-begin[3])/60) 
    }