rdirichlet_function(param) {
  tmp_rgamma(length(param), param)
  tmp/sum(tmp) }

##############################################################################
### MCMC Sampler for two-pop, two-test analysis, without assuming independence
### of tests, i.e. DUAL-IND of Ch. 5.
### LAST UPDATE: Mar. 10/04
###
### Beta priors on sens. (p1 and p2) and spec. (q1 and q2),
### hyperparamters a.XX and b.XX
###
### conditional, truncated exponential priors on covariance between tests
### for truly exposed (dp) and truly unexposed (dq), as described in Ch. 5,
### k is the hyperparameter governing these priors
### 
### Uniform priors on prevalences (ra and rb).
###
### Input is two 2*2 matrices, exposure assessments for controls and cases,
### i.e. row describes first assessment (unexposed/exposed), col describes 2nd
###
### Output is keep * 9 matrix, cols are sampled vals of 
### (p1,p2,dp,q1,q2,dq,ra,rb,log-OR) 
##############################################################################


dual.dep <- function(tab.con, tab.cas, 
              a.p1=1, a.p2=1, a.q1=1, a.q2=1, 
              b.p1=1, b.p2=1, b.q1=1, b.q2=1,
              k=4*log(4),
              burn=1000,keep=5000) {

### initial values
p1 <- p2 <- q1 <- q2 <- 0.9; dp <- 0.01; dq <- 0.01
r1 <- tab.con[2,2]/(tab.con[1,1]+tab.con[2,2])
r2 <- tab.cas[2,2]/(tab.cas[1,1]+tab.cas[2,2])
tab.xps <- tab.unx <- matrix(NA,2,2)
res <- matrix(0,keep,9)
dimnames(res) <- list(NULL, c("p1","p2","dp","q1","q2","dq","r1","r2","psi"))

for (mnlp in ((-burn):keep)) {

  if ((mnlp%%100)==0) {print(mnlp)}

  ### update exps
  cnt1 <- cnt2 <- 0

  ### (1,1) --   
  pr.x <- (1-p1)*(1-p2)+dp; pr.u <- q1*q2+dq
  tmp1 <- rbinom(1,size=tab.con[1,1], prob=pr.x*r1/(pr.x*r1+pr.u*(1-r1)) )
  tmp2 <- rbinom(1,size=tab.cas[1,1], prob=pr.x*r2/(pr.x*r2+pr.u*(1-r2)) )
  cnt1 <- cnt1+tmp1; cnt2 <- cnt2+tmp2
  tab.xps[1,1] <- tmp1+tmp2
  tab.unx[1,1] <- tab.con[1,1]+tab.cas[1,1]-tab.xps[1,1]

  ### (1,2) -+   
  pr.x <- (1-p1)*p2-dp; pr.u <- q1*(1-q2)-dq
  tmp1 <- rbinom(1,size=tab.con[1,2], prob=pr.x*r1/(pr.x*r1+pr.u*(1-r1)) )
  tmp2 <- rbinom(1,size=tab.cas[1,2], prob=pr.x*r2/(pr.x*r2+pr.u*(1-r2)) )
  cnt1 <- cnt1+tmp1; cnt2 <- cnt2+tmp2
  tab.xps[1,2] <- tmp1+tmp2
  tab.unx[1,2] <- tab.con[1,2]+tab.cas[1,2]-tab.xps[1,2]  

  ### (2,1) +-   
  pr.x <- p1*(1-p2)-dp; pr.u <- (1-q1)*q2-dq
  tmp1 <- rbinom(1,size=tab.con[2,1], prob=pr.x*r1/(pr.x*r1+pr.u*(1-r1)) )
  tmp2 <- rbinom(1,size=tab.cas[2,1], prob=pr.x*r2/(pr.x*r2+pr.u*(1-r2)) )
  cnt1 <- cnt1+tmp1; cnt2 <- cnt2+tmp2
  tab.xps[2,1] <- tmp1+tmp2
  tab.unx[2,1] <- tab.con[2,1]+tab.cas[2,1]-tab.xps[2,1]  

  ### (2,2) ++   
  pr.x <- p1*p2+dp; pr.u <- (1-q1)*(1-q2)+dq
  tmp1 <- rbinom(1,size=tab.con[2,2], prob=pr.x*r1/(pr.x*r1+pr.u*(1-r1)) )
  tmp2 <- rbinom(1,size=tab.cas[2,2], prob=pr.x*r2/(pr.x*r2+pr.u*(1-r2)) )
  cnt1 <- cnt1+tmp1; cnt2 <- cnt2+tmp2
  tab.xps[2,2] <- tmp1+tmp2
  tab.unx[2,2] <- tab.con[2,2]+tab.cas[2,2]-tab.xps[2,2]      

  ### update r1, r2
  r1 <- rbeta(1, cnt1+1, sum(tab.con)-cnt1+1)
  r2 <- rbeta(1, cnt2+1, sum(tab.cas)-cnt2+1)

  ### update p1,p2,dp
  tht <- c(p1*p2+dp, p1*(1-p2)-dp, (1-p1)*p2-dp)
  tht.str <- rdirichlet(c(1,1,1,1)+ 
    c(tab.xps[2,2],tab.xps[2,1],tab.xps[1,2],tab.xps[1,1]))[1:3]
  p1.str <- tht.str[1]+tht.str[2]; p2.str <- tht.str[1]+tht.str[3]
  dp.str <- tht.str[1]*(1-sum(tht.str))-tht.str[2]*tht.str[3]

  if ( (dp.str>0) && ((p1.str+q1)>1) && ((p2.str+q2)>1) ) {
    invscl <- k/sqrt(p1.str*(1-p1.str)*p2.str*(1-p2.str))
    maxval <- min(p1.str,p2.str)-p1.str*p2.str
    acc <- invscl*exp(-invscl*dp.str)/(1-exp(-invscl*maxval))
    invscl <- k/sqrt(p1*(1-p1)*p2*(1-p2))
    maxval <- min(p1,p2)-p1*p2
    acc <- acc*(1-exp(-invscl*maxval)) / (invscl*exp(-invscl*dp)) 
    if (runif(1)<acc) {
      p1 <- p1.str; p2 <- p2.str; dp <- dp.str } }

  ### update q1,q2,dq
  tht <- c(q1*q2+dq, q1*(1-q2)-dq, (1-q1)*q2-dq)
  tht.str <- rdirichlet(c(1,1,1,1)+
    c(tab.unx[1,1], tab.unx[1,2], tab.unx[2,1], tab.unx[2,2]))[1:3]
  q1.str <- tht.str[1]+tht.str[2]; q2.str <- tht.str[1]+tht.str[3]
  dq.str <- tht.str[1]*(1-sum(tht.str))-tht.str[2]*tht.str[3]
  if ( (dq.str>0) && ((p1+q1.str)>1) && ((p2+q2.str)>1) ) {
    invscl <- k/sqrt(q1.str*(1-q1.str)*q2.str*(1-q2.str))
    maxval <- min(q1.str,q2.str)-q1.str*q2.str
    acc <- invscl*exp(-invscl*dq.str)/(1-exp(-invscl*maxval))
    invscl <- k/sqrt(q1*(1-q1)*q2*(1-q2))
    maxval <- min(q1,q2)-q1*q2
    acc <- acc*(1-exp(-invscl*maxval)) / (invscl*exp(-invscl*dq)) 
    if (runif(1)<acc) {
      q1 <- q1.str; q2 <- q2.str; dq <- dq.str } }

  if (mnlp>0) {
    res[mnlp,] <- c(p1,p2,dp,q1,q2,dq,r1,r2,log(r2/(1-r2))-log(r1/(1-r1))) } }

res }