#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))