#################################################### ### Gibbs Sampler for two-pop, two-test analysis ### ### Beta priors on p1, p2, q1, q2 ### ### Unif prior on ra, rb ### #################################################### fitB_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, burn=1000,keep=5000) { ### initial values p1_p2_q1_q2_0.95 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,7) dimnames(res)_list(NULL, c("p1","p2","q1","q2","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); pr.u_q1*q2 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; pr.u_q1*(1-q2) 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); pr.u_(1-q1)*q2 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; pr.u_(1-q1)*(1-q2) 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 tmp_as.vector(tab.xps[,1]+tab.xps[,2]) flg_T; while (flg) { p1_rbeta(1, tmp[2]+a.p1, tmp[1]+b.p1); flg_((p1+q1)<1) } ### update p2 tmp_as.vector(tab.xps[1,]+tab.xps[2,]) flg_T; while (flg) { p2_rbeta(1, tmp[2]+a.p2, tmp[1]+b.p2); flg_((p2+q2)<1) } ### update q1 tmp_as.vector(tab.unx[,1]+tab.unx[,2]) flg_T; while (flg) { q1_rbeta(1, tmp[1]+a.q1, tmp[2]+b.q1); flg_((p1+q1)<1) } ### update q2 tmp_as.vector(tab.unx[1,]+tab.unx[2,]) flg_T; while (flg) { q2_rbeta(1, tmp[1]+a.q2, tmp[2]+b.q2); flg_((p2+q2)<1) } if (mnlp>0) { res[mnlp,]_c(p1,p2,q1,q2,r1,r2,log(r2/(1-r2))-log(r1/(1-r1))) } } res }