### PooledHIVTest

################################################################################# 

### This program calculates the probability that a single blood sample contains
### HIV given that it was pooled with other blood samples and the pooled sample
### tests positive.

### The program specifies n = the number of blood samples that are in a single
### pool and reps = the total number of (individual) blood samples that will be
### simulated.  For example, if n = 5 and reps = 100 then there are 20 pools of 
### size 5 simulated.

### In the function PoolTest a pool of size n is simulated and the NumWtihHIV is 
### recorded.  Then if at least one (simulated) person has HIV then the test has
### a probability of 0.98 of detecting HIV.  If no person has HIV then the test
### has probability 0.005 of mistakenly dtecting HIV.  The counters Positives and 
### TruePostives are incremented.  If the test is positive the Positives are
### incremented by n and TruePositive are incremented by NumWith HIV.

### Lastly the simulated conditional probability is reported:
###	 p = TruePostives / Postives

### The actual proabilities work out to be:
### 	n = 5,  p = 0.1935
### 	n = 10, p = 0.1036
### 	n = 20, p = 0.0567

#################################################################################

# function that simulates poolled testing.

TestPool <- function(n,Positives,TruePositives){
	
	# Simulate the blood samples.
	
	NumWithHIV <- 0
	for(j in 1:n){
		if(runif(1) < 0.015) NumWithHIV <- NumWithHIV + 1
	}
	#print(NumWithHIV)

	# Simulate the blood test of the pooled sample.

	TestPositive <- 0
	if(NumWithHIV >= 1){
		if(runif(1) <= 0.98) TestPositive <- 1
	} else{
		if(runif(1) <= 0.005) TestPositive <- 1
	}
	#print(TestPositive)

	# Update counters.
	
	if(TestPositive > 0){
		Positives <- Positives + n
		TruePositives <- TruePositives + NumWithHIV
	}
	#print(Positives)
	#print(TruePositives)
	
	# Return results.
	
	result <- c(Positives,TruePositives)
	return(result)
}

# pool size.

n <- 20

# a.

# total number of people to be tested.  this number should be a multiple of n.

reps <- 100

Positives <- 0
TruePositives <- 0

i <- 0
repeat{
	result <- TestPool(n,Positives,TruePositives)
	Positives <- result[1]
	TruePositives <- result[2]

	i <- i + n
	if(i >= reps) break
}

Positives <- result[1]
TruePositives <- result[2]

p <- TruePositives/Positives

cat("For pool size",n,"and",reps,"people\n")
cat("the approximate proability is",p,".\n")

# b.
        
# total number of people to be tested.  this number should be a multiple of n.
        
reps <- 1000

Positives <- 0
TruePositives <- 0
        
i <- 0
repeat{  
        result <- TestPool(n,Positives,TruePositives)
        Positives <- result[1]
        TruePositives <- result[2]
        
        i <- i + n
        if(i >= reps) break
}
         
Positives <- result[1]  
TruePositives <- result[2]  
        
p <- TruePositives/Positives

cat("For pool size",n,"and",reps,"people\n")
cat("the approximate proability is",p,".\n")

# c.

# total number of people to be tested.  this number should be a multiple of n.

reps <- 10000

Positives <- 0
TruePositives <- 0

i <- 0
repeat{
        result <- TestPool(n,Positives,TruePositives)
        Positives <- result[1]
        TruePositives <- result[2]

        i <- i + n
        if(i >= reps) break
}       
        
Positives <- result[1]
TruePositives <- result[2]
        
p <- TruePositives/Positives
 
cat("For pool size",n,"and",reps,"people\n")
cat("the approximate proability is",p,".\n")