Chapter 2 Hello bookdown

All chapters start with a first-level heading followed by your chapter title, like the line above. There should be only one first-level heading (#) per .Rmd file.

2.1 A section

All chapter sections start with a second-level (##) or higher heading followed by your section title, like the sections above and below here. You can have as many as you want within a chapter.

An unnumbered section

Chapters and sections are numbered by default. To un-number a heading, add a {.unnumbered} or the shorter {-} at the end of the heading, like in this section.

############################## PART A ########################################

# We need to use Monte Carlo integration here. In other words, we'll do
# these steps: 
# 1. get a vector of expected losses for a loss function 
# 2. Get a bunch more vectors. Like 1000. 
# 3. Take the average of all those expected loss vectors 
# 4. Plot that average vs your actions. 

n <- 100 # sample size
mc <- 100 # number of Monte Carlo runs

# calculate gamma density for each x value
alpha <- 15  # shape parameter
beta  <- 3   # rate parameter

# y is sampled from a gamma distribution with parameters alpha and beta.
# Gives us a 1 x n vector of y randomly sampled y values from this
# distribution.
yRandVec <- rgamma(n, alpha, beta)

# create a vector of actions. We'll eventually calculate a distinct
# expected loss for each one of these actions:
a = 10
actions <- seq(a/n, a, by=a/n )
# cut off the last element to ensure actions has 1000 elements:
#actions <- actions[1:(length(actions)-1)]

# Loop through and calculate the loss for every possible action. After 
# the loop, we end up with a 1 x n vector of expected losses.
emptyValuesForMatrix <- rep(NA, mc*n)
expectedLossesMatrix <- matrix(emptyValuesForMatrix, nrow = mc, ncol = n)

for (iter0 in 1:mc)
{
  expectedLosses <- vector(mode = "double", length = n)
  for (iter in 1:n)
  {
    expectedLosses[iter] = mean( (actions[iter] - yRandVec)^2 )
  }
  # Now store that expected loss vector in our giant matrix: 
  expectedLossesMatrix[iter0, ] = expectedLosses
}

mcExpectedLosses = colMeans(expectedLossesMatrix)
plot(actions, mcExpectedLosses)


# Find the minimizer of the loss function-- that is, what action gives us 
# the smallest value of expectedLosses?
indexOfSmallestLoss <-  which.min(mcExpectedLosses)
minimizerAction <-  actions[indexOfSmallestLoss]
cat("Minimizer of expected loss function: ", minimizerAction, "\n")


# We've now confirmed that the loss minimizer is equal to the mean of Y, or
# alpha/beta. Now, we'll confirm that the value of the expected loss itself at 
# at the minimizer is the variance of Y, or alpha/(beta^2).
smallestLoss <-  mcExpectedLosses[[indexOfSmallestLoss]] 
# turning into a list for some reason
cat("Value of expected loss at optimal action: ", smallestLoss, "\n")


################################ PART B ######################################

# We'll now computed the expected loss for the check loss function. 
# The check loss function is computed differently depending on the value
# of y, so we need to include logic to perform the correct calculation.

expectedCheckLosses <- vector(mode = "double", length = n)
checkLosses <- vector(mode = "double", length = length(yRandVec))
emptyValuesForMatrix <- rep(NA, mc*n)
expectedCheckLossesMatrix <- matrix(emptyValuesForMatrix, nrow = mc, ncol = n)
q = 0.5

for (iter0 in 1:mc)
{
  for (iter in 1:n)
  {
    # Check y vs a for every value of y
    for (j in 1:length(yRandVec))
    {
      # If the realization of the random variable Y is greater than a, compute
      # this form of the check loss function (at this j):
      if ((yRandVec[j] > actions[iter]) == TRUE)
      {
        checkLosses[j] = q*(yRandVec[j] - actions[iter])  
        
      }
      # Otherwise, use the other form of the check loss function (at this j):
      else
      {
        checkLosses[j] = (1-q)*(actions[iter] - yRandVec[j])
      }
    }
    expectedCheckLosses[iter] = mean(checkLosses)
  }
  # Now fill in our matrix of expected losses, just like last time:
  expectedCheckLossesMatrix[iter0, ] = expectedCheckLosses
}

mcExpectedCheckLosses = colMeans(expectedCheckLossesMatrix)

plot(actions, mcExpectedCheckLosses)
Finvq = qgamma(q, alpha, beta)

cat("\n")
cat("Value of q: ", q, "\n")
cat("Optimal action of expected check loss function",
    actions[which.min(mcExpectedCheckLosses)], "\n" )
cat("Value of Gamma inverse cdf at q: ", Finvq, "\n")


############################# PART C #########################################

tau <- 5.3
yRandVecPoisson <- rpoisson(n, tau)

# We're using an indicator function. That is, indicator(y) = 1 if 
# y == a and 0 otherwise. 

for (iter in 1:n)
{
  
}