I am trying to simulate the answer to this problem in R.
Northwest Minnesota is prone to flooding in the Spring. Suppose that 20% of all farmers are insured against flood damage. Four farmers are selected at random. What is the probability that at least two farmers have flood insurance?
Thanks for any help from anyone.
This answer walks you through steps that are common to most stochastic simulations, showing how to create the code in small, simple, easily-tested chunks on any software platform. The process is illustrated with R code. You can run the code snippets as you go along to see what they produce.
FWIW, here's a compact, quick-and-dirty solution that compresses the calculations into one line. It prints out an estimate based on a simulation with n iterations. (Figure around four million iterations per second of CPU time.)
n <- 1e6; mean(rowSums(matrix(runif(4*n) <= 0.20, nrow=n)) >= 2)
Create a box with tickets, one per farmer. On those tickets are written "insured" and "not insured." The proportion of "insured" in the box must be 20%. (For a detailed account of the tickets-in-a-box model of random variables, see What is meant by a random variable
In almost every software system that supports random number generation, there is a function to draw a ticket out of a very large box (and replace it after the ticket is observed). These tickets have floating point values from $0$ to just a tiny bit less than $1$ written on them. For any interval $0 \le a \le b \lt 1$, the proportion of tickets with values between $a$ and $b$ is $b-a$. You can exploit this function by writing "insured" on 20% of all the numbers between $0$ and $1$. For instance, you could write "insured" on all numbers between $0$ and $20/100$.
This is the Uniform distribution (between $0$ and $1$). In R this function is called runif. A single draw of a ticket from this box could be programmed as
label <- ifelse(runif(1) <= 20/100, "insured", "not insured")
It is faster and more convenient to dispense with the label, though, and replace it with the number $1$ (reserving $0$ for all other tickets). This gives the simpler R code
label.indicator <- runif(1) <= 20/100
because in R (as in many systems) a true value is equated with $1$ and a false value with $0$.
Just draw tickets from the box independently. That's done with a loop. In R the loop is performed for you (very efficiently, behind the scenes) when you request multiple values from runif. What is important to know is that these values are (pseudo) independent: they do not appear to depend on one another. Thus,
label.indicators <- runif(4) <= 20/100
will simulate drawing the tickets of four farmers from this box. It produces an array of four numbers from the set {FALSE, TRUE} or equivalently $\{0,1\}$.
The number of farmers in any group is found by summing their indicators.
farmer.count <- sum(label.indicators)
This is because each insured farmer contributes a $1$ to the sum and each uninsured farmer contributes a $0$. The sum merely counts the insured farmers. The efficiency of using indicators, instead of labels, is apparent here.
This is a loop in most platforms. In some, including R, it's faster to draw lots of tickets and group them in sets of four. (This is because often it's almost as quick to generate many random values as it is to generate just one: there's less overhead involved.) Each group represents one iteration of the simulation. The following code puts the four tickets for each iteration in rows of an array and then sums each row (as in Step 3).
n.sim <- 1e6 # Number of iterations
n.farmers <- 4 # Number of farmers per iteration
simulation <- rowSums(matrix(runif(n.sim * n.farmers) <= 20/100, nrow=n.sim))
The chance of observing two or more insured farmers can be estimated as the proportion of iterations in which two or more insured farmers were found in the sample. As before, this is efficiently found by testing and summing (or testing and averaging):
estimate <- mean(simulation >= 2)
You have just performed a computer experiment. Like any other experiment, the results are variable, so you ought to provide a standard error for the result. In this Binomial experiment, the standard error of the estimate $\hat p$ (with $n$ iterations) is
$$\text{se}(\hat p) = \sqrt{\hat p\left(1 - \hat p\right) / n}.$$
Compute this and print out the results:
se.estimate <- sqrt(estimate * (1-estimate) / n.sim)
print(c(Estimate=estimate, SE=se.estimate), digits=2)
When the random seed is set to $17$ (see the full code below) and a million iterations are run, the output is
Estimate SE
0.18127 0.00039
Here's the full solution. It is structured to permit easy variation of the parameters so, by repeating it (the calculation takes less than a second), you can study how the results depend on the parameters to get a more intuitive feel for what is happening.
#
# Specify the problem.
#
p <- 20/100 # Chance of insurance
k <- 2 # Minimum number of insured farmers
n.farmers <- 4 # Number of farmers per iteration
n.sim <- 1e6 # Number of iterations
#
# Simulate.
#
set.seed(17) # Optional: provides a reproducible result
simulation <- rowSums(matrix(runif(n.sim * n.farmers) <= p, nrow=n.sim))
#
# Post-process and report.
#
estimate <- mean(simulation >= k)
se.estimate <- sqrt(estimate * (1-estimate) / n.sim)
print(c(Estimate=estimate, SE=se.estimate), digits=2)
Here is a possible solution to your question
prop=0
T=1e6
for (t in 1:T){
insrd=sample(0:1,4,rep=TRUE,prob=c(8,2))
prop=prop+(sum(insrd)>1)}
print(prop/T)
with answer 0.1809. But you can compute the exact value [0.1808] of this probability by realising that the draw of four farmers is a Binomial B(4,0.2) random variable when counting the number of insured farmers out of the four.
