# Example simulation program for Exercise 2.5-3 (Example 2.5-4) from
# the Higgins & Keller-McNulty text.


# This function simulates n values for the random variable W (the sum of lost 
# revenues for five flights), and returns the empirical probability mass 
# function as a list of two vectors, giving the possible values and the 
# estimated probability for each value.  Note that the possible values are
# 0, 70, 140, etc. up to 5*210=1050.
# 
# This function isn't written in the most efficient way.  It's meant to
# illustrate the idea.

simulate.flights <- function (n)
{
  # Find the set of possible values for W, which form a sequence from 0
  # to 1050 in steps of 70.

  possible.values <- seq(0,1050,by=70)

  # Initialize the counts (as many as there are possible values) to zero.

  counts <- rep(0,length(possible.values))

  # Main simulation loop.  Simulates one week of five flights each iteration.

  for (i in 1:n)
  { 
    # Simulate a value for W, by adding together five simulated values for Y.

    w <- 0
    for (j in 1:5)
    { u <- runif(1)
      if (u<0.5) y <- 0
      else if (u<0.8) y <- 70
      else if (u<0.95) y <- 140
      else y <- 210
      w <- w+y
    }

    # Find where this value for w is located in the w and p vectors.

    k <- 1
    while (possible.values[k]!=w)
    { k <- k+1
      if (k>length(possible.values))
      { stop("Got an impossible value for w!")
      }
    }

    # Increment the count of how many times this value has occurred.

    counts[k] <- counts[k] + 1
  }

  # Divide the counts by number of flights simulated to get the 
  # estimated probabilities.  This is a vector operation, in which
  # every element of "count" is divided by n.

  probabilities <- counts / n

  # Return the possible values and their estimated probabilities.

  list (possible.values=possible.values, probabilities=probabilities)
}