Mostly we have been estimating probabilities.

• We repeated the chance process a very large number of times
• On each repetition, we recorded whether or not the event occurred.
• After we finished all of the repetitions, we computed the proportion of those repetitions in which the event did occur.
• This proportion is our estimate of the probability of the event.
• The more repetitions we arrange for, the closer our estimate is liable to be to the actual probability of the event.

Let's switch gears: estimate the expected value of a random variable.

Recall:

• a random variable is a number whose value is the result of some chance process.
• the expected value of the random variable is the average of its values over many, many repetitions of that chance process.

Consider this chance process:

• Pick a real number \(x_1\) between 0 and 1 at random.
• Pick another real number \(x_2\) at random (also between 0 and 1). Add it to \(x_1\).
• If the sum \(x_1 + x_2\) exceeds 1, then stop.
• Otherwise, pick a third real number \(x_3\), and add it to the sum of the previous two.
• If the sum \(x_1+x_2+x_3\) exceeds 1, then stop.
• Otherwise, pick a fourth real number …
• Keep going until the sum of the numbers picked exceeds 1.

Let \(X\) be the number of numbers you have to pick.

Here is a function to simulate \(X\):

numberNeeded <- function(verbose = FALSE) {
  mySum <- 0
  count <- 0
  while( mySum < 1 ) {
    number <- runif(1)
    mySum <- mySum + number
    count <- count + 1
    if ( verbose ) {
      cat("Just now picked ", number, ".\n", sep = "")
      cat(count, " number(s) picked so far.\n", sep = "")
      cat("Their sum is ", mySum, ".\n\n", sep = "")
    }
  }
  count
}

numberNeeded(verbose = TRUE)

Try it a few times.

needed <- numeric(1000)
for (i in 1:1000 ) {
  needed[i] <- numberNeeded()
}
cat("The expected number needed is about ", 
    mean(needed), ".\n", sep = "")

## The expected number needed is about 2.697.

table(needed)

## needed
##   2   3   4   5   6   7 
## 528 303 123  37   8   1

prop.table(table(needed))

## needed
##     2     3     4     5     6     7 
## 0.528 0.303 0.123 0.037 0.008 0.001

This estimates the distribution of \(X\), the number needed.

First, generalize numberNeeded() to handle target-numbers other than 1:

numberNeeded <- function(target) {
    mySum <- 0
    count <- 0
    while( mySum < target ) {
      number <- runif(1)
     mySum <- mySum + number
      count <- count + 1
    }
    count
  }

numberNeededSim <- function(target = 1, reps = 1000, 
                            seed = NULL, report = FALSE) {
  if ( !is.null(seed) ) {
    set.seed(seed)
  }
  needed <- numeric(reps)   # set up vector for results
  for (i in 1:reps ) {      # loop to repeat random process
    needed[i] <- numberNeeded(target)   # store result each time through
  }
  if ( report ) {
    print(prop.table(table(needed)))
    cat("\n")
    cat("The expected number needed is about ",
        mean(needed), ".\n", sep = "")
  }
  invisible(mean(needed))
}

numberNeededSim(target = 1, reps = 10000, seed = 4848, report = TRUE)

## needed
##      2      3      4      5      6      7      8 
## 0.4937 0.3350 0.1313 0.0317 0.0071 0.0010 0.0002 
## 
## The expected number needed is about 2.7273.

It is known that when the target number is 1 the expected number of numbers needed is:

\[\mathrm{e} \approx 2.71828 \ldots\]

Here is the template for a function to estimate a probability via Monte Carlo Simulation:

probEstimation <- function(reps) {
  results <- logical(reps)  # results will go here
  for ( 1 in 1:reps ) {
    # simulate the random process,
    # getting outcome (TRUE if event happens, FALSE if not)
    results[i] <- outcome  # store the outcome
  }
  mean(results)  # return estimate of probability of event
}

Here is the template for a function to estimate an expected value via Monte Carlo Simulation:

expValEstimation <- function(reps) {
  results <- numerical(reps)  # results will go here
  for ( 1 in 1:reps ) {
    # simulate the random process,
    # getting a value for the random variable
    results[i] <- value  # store the value
  }
  mean(results)  # return estimate of the expected value
}

Let's cover this from the textbook.

https://homerhanumat.github.io/r-notes/appeals-court-paradox.html

https://homerhanumat.github.io/r-notes/drunken-turtle-sim.html