#STATS: Gibbs Sampling, Sec 8. #C. Fraser 7/10/99 #Gibbs Sampler Program set.seed(1) #Set initial value for random number #generator to obtain reproducible #results. Omit for varying random #simulation runs. #Known values: Sensitivity (eta) and specificity (theta) eta <- .99 #Sensitivity theta <- .97 #Specificity #Prior density parameters alpha <- 2 #Prevalence prior parameter beta <- 50 #Prevalence prior parameter #Initial test values M2 <- 20000 #Number of iterations M1 <- M2/2 #Number of burn-in iterations A <- 49 #DATA: Number testing postive B <- 1000 - A #DATA: Number testing negative X <- numeric(M2) #Number of true positives Y <- numeric(M2) #Number of false positives pi <- numeric(M2) #Prevalence given A, B, X, Y Sum <- 0 #Running sum of sampled values of pi Average <- numeric(M2) #Running averages pi[1] <- 0 #Initialize prevalence vector Average[1] <- pi[1] #Initialize average vector #[Vectors start at m=1] #Generate the chain of prevalence values for (m in 2:M2) { X[m] <- rbinom(1, A, (pi[m-1]*eta)/ ((pi[m-1]*eta)+(1-pi[m-1])*(1-theta))) Y[m] <- rbinom(1, B, (pi[m-1]*(1-eta))/ ((pi[m-1]*(1-eta))+(1-pi[m-1])*theta)) pi[m] <- rbeta(1, (X[m] + Y[m] + alpha), (A + B - X[m] - Y[m] + beta)) Sum <- Sum + pi[m] Average[m] <- Sum / m } #Compute quantiles of histogram Gibbs.interval <- quantile(pi[(M1+1):M2], c(.025, .975)) #Compute mean of pi values after burn-in mu <- mean(pi[(M1+1):M2]) #Compute median of pi values after burn-in median <- median(pi[(M1+1):M2]) #Plot histogram of sampled values of pi (after burn-in) hist(pi[(M1+1):M2], main = "Figure 3. Sampled Prevalence Values (After Burn-In).", xlab = "Prevalence", ylab = "Frequency", ylim = c(0, 3000)) #Set page format #par(mfcol=c(1,2)) #Plot all sampled values of pi vs. iteration plot(pi, main = "Figure 4. Sampled Prevalence Values in Sequence.", xlab = "Iteration", ylab = "Prevalence", type = "l", ylim = c(0, .1)) #Plot all running averages of pi vs. iteration plot(Average[2:M2], main = "Figure 5. Running Averages of Sampled Prevalences.", xlab = "Iteration", ylab = "Average", type = "l", ylim = c(.015, .025)) #Print results Gibbs.interval mu median