#STATS: Gibbs Sampling, Sec 4.
#C. Fraser 7/10/99

#Gibbs Matrix Program

set.seed(1)                         #Set initial value for random number 
                                    #generator to obtain reproducible 
                                    #results. Omit for varying random 
                                    #simulation runs.

M2 <- 50000                         #Total number of iterations to be performed 
M1 <- M2 / 2                        #Number of iterations in burn-in period 

eta     <- .99                      #Sensitivity of diagnostic test
theta   <- .97                      #Specificity of diagnostic test 
gamma   <- .4024                    #Predictive value of a positive test
delta   <- .9998                    #Predictive value of a negative test


Disease <- numeric(M2)              #Vector of disease values (0 or 1)
Test    <- numeric(M2)              #Vector of test values (0 or 1)
Pi.prop <- numeric(M2)              #Vector of proportions (testing positive)
Sum     <- 0                        #Running sum (initialized to zero)
Disease[1] <- 0                     #Initialize disease vector
                                    #[Starts at m=1]
            
#Generate Values
for (m in 2:M2)
{      
      if (Disease[m-1] == 1)
            Test[m-1] <- rbinom(1,1,eta)
      else
            Test[m-1] <- rbinom(1,1,1-theta)
                                                                                    
      if (Test[m-1] == 1) 
            Disease[m] <- rbinom(1,1,gamma)
      else
            Disease[m] <- rbinom(1,1,1-delta)
                  
      Sum <- Sum + Disease[m]
      Pi.prop[m] <- Sum / m
}

#Calculate prevalence estimate
D1count <- sum(Disease[(M1+1):M2])            
Pi.post <- D1count / (M2 - M1)

#Output estimated prevalence
Pi.post

#Plot of running proportion of infected samples
plot(Pi.prop, main = "Figure 1. Running Proportion of Infected Samples.",
     xlab = "Iteration", ylab = "Proportion", type = "l", ylim = c(0, .03))