#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