The following solves a decision problem taken from BDA3 Ch9 exercise 1, with slight modifications.

The production decision problem

The goal is to maximize expected net profits from producing widgets. The widgets cost $2 to make and sell for $3. A forecast for the market for widgets provides a source of uncertainty. The forecast prediction is normally distributed with mean of 10,000 and sd of 5,000. So how many widgets should you make?

I’m going to follow the basic steps outlined by Gelman et al. (2013), but there are probably faster ways to calculate the optimal decision.

First, define the set of decisions. Because the decision determines production this will also denote the set of outcomes (x) which implies there’s no uncertainty yet.

d <- seq(from = 1, to = 20000, by = 1)

Next, we will define the utility function, which in this example is just a net revenue function. If you produce more widgets than the market demands these cannot be sold, but you still pay the cost of production. I’ll code this as function which will be useful later.

u <- function(x, m) {
  # This function takes two arguments: the decision of how much to produce (x), and the market
  # size (m).
  profit <- 3 * min(x, m) - 2 * x
  return(profit)
}

Now lets calculate the expected utility of producing 1000 widgets given 10000 draws from our forecasted market. Remember, these forecasters are clever folks, and truncated the normal distribution at 0 by throwing out values less than that.

set.seed(1234)

m_for <- rnorm(10000, mean = 10000, sd = 5000)
m_for <- m_for[m_for > 0]

hist(m_for)

Given these draws the expected profit is…

x <- 1000 # Our decision guess.

exp_profit <- function(x, m_for) {
  profit <- vector(length =  length(m_for))
  for(i in 1:length(m_for)) {
    profit[i] <- u(x, m_for[i])
  }
#  hist(profit)
  ex_prof <- sum(profit)/length(m_for)
  return(ex_prof)
}

exp_profit(1000, m_for)
## [1] 981.8439

My priors tell me we can do better. Let’s optimize how many widgets to produce.

max_d <- c(0,0) # An "empty" vector indicating optimal decision and expected profit.

for(i in 1:length(d)) {
  p <- exp_profit(d[i], m_for)
  if(p > max_d[2]) 
    max_d <- c(d[i], p)
  else max_d <- max_d
}

max_d
## [1] 8150.000 5190.505

And there we have it!

So, how much would you pay those clever forecasters for a better forecast? Let’s imagine they can reduce the standard deviation to 2,500. What would that information be worth?

Here’s their new and improved forecast.

m_for2 <- rnorm(10000, mean = 10000, sd = 2500)
m_for2 <- m_for2[m_for2 > 0]

Let’s solve the optimization.

max_d2 <- c(0,0) # An "empty" vector indicating optimal decision and expected profit.

for(i in 1:length(d)) {
  p <- exp_profit(d[i], m_for2)
  if(p > max_d2[2]) 
    max_d2 <- c(d[i], p)
  else max_d2 <- max_d2
}

max_d2
## [1] 8879.000 7224.588

The improved information increases your profits by $2,034.08, so the new forecast is worth at most that much.