In cboettig/multiple_uncertainty: Multiple Uncertainty

library("stringr")
library("dplyr")
library("tidyr")
library("ggplot2")
library("lazyeval")
library("readr")
library("multipleuncertainty")
knitr::opts_chunk$set(cache = TRUE, comment=NA)

Optimal control under multiple uncertainty by model:

mean_value <- function(sigma_g_belief, sigma_m_belief, sigma_i_belief,
                     sigma_g_true, sigma_m_true, sigma_i_true, f, noise_dist){

  delta <- 0.05
  grid <- seq(0,150,length=151) #seq(0, 200, length = 401)

  S <- multiple_uncertainty(f = f, 
                            x_grid = grid, 
                            sigma_g = sigma_g_belief, 
                            sigma_m = sigma_m_belief, 
                            sigma_i = sigma_i_belief,
                            delta = delta,
                            noise_dist = noise_dist)

  ## Function to compute to do a replicate simulation and return NPV
  replicates <- function(){
    sim <- scenario(S, f = f, 
                    x_grid = grid, 
                    sigma_g = sigma_g_true, 
                    sigma_m = sigma_m_true, 
                    sigma_i = sigma_i_true, 
                    noise_dist = noise_dist)
    sim %>% 
      mutate(value = h * (1-delta) ^ t) %>% 
      summarise(NPV = sum(value))
  }

  ## Simulate 100 replicates of the policy and return E(NPV)
  data.frame(rep = 1:100) %>% 
    group_by(rep) %>% 
    do(replicates()) %>% 
    ungroup() %>%
    summarise(ENPV = mean(NPV), sd = sd(NPV))


  }

low <- 0.01
high <- 0.5
cases <-
  expand.grid(
  sigma_g_belief = c(low, high),
  sigma_m_belief = c(low, high),
  sigma_i_belief = c(low, high),
  sigma_g_true = c(low, high),
  sigma_m_true = c(low, high),
  sigma_i_true = c(low, high)
  )

group_by_all <- function(x,...) x %>%  group_by_(.dots = lapply(names(x), as.name))  
cases %>%
  group_by_all() %>%
  do(mean_value(.$sigma_g_belief, .$sigma_m_belief, .$sigma_i_belief,
                 .$sigma_g_true, .$sigma_m_true, .$sigma_i_true,
                 "ricker", "lognormal")) %>% 
  ungroup() -> df

We want to collapse the columns that simply list the values of each true/belief sigma value into a single scenario key, where we'll use L for "Low" noise and H for high. We indicate the belief values of growth, measurement, and implementation error, respectively, followed by the true values used for the simulations. Thus HLL | LHL indicates that the optimal solution was computing assuming high growth error, but low measurement and implementation, and then this policy was applied to simulations in which growth noise and implementation noise were low, but measurement was high.

rekey_and_normalize <- function(df, low, high){

  ## 
  rekey <- function(str){
    str %>% 
      stringr::str_replace_all(as.character(low), "L") %>%
      stringr::str_replace_all(as.character(high), "H") %>%
      stringr::str_replace_all("_", "")
  }

  ## Replace columns with noise values like (0.5, 0.5, .1) into keys like (HHL)
  df %>% 
    unite_("true", names(df)[4:6]) %>% 
    mutate(true = rekey(true)) %>% 
    unite_("belief", names(df)[1:3]) %>% 
    mutate(belief = rekey(belief)) %>% 
    select(belief, true, ENPV, sd) -> 
    table

  ## Include unique id for binding
  table <- cbind(id = 1:dim(df)[1], table)

  ## Normalize
  table %>% filter(belief == true) -> optimal
  tbls <- lapply(optimal$true, function(x){
    expr1 <- interp(~true == x)
    expr2 <- interp(~N, N = filter_(optimal, .dots = list(expr1))[["ENPV"]])
    table %>% filter_(~true == x) %>% mutate_(.dots = setNames(list(expr2), "N"))
  })
  tbl <- bind_rows(tbls) %>% arrange(id)

  tbl
}

tbl <- rekey_and_normalize(df, low, high)

In these 8 scenarios, the beliefs match the true values for each of the noise levels:

tbl %>% filter(belief == true) %>% knitr::kable("pandoc")

Summary and analysis

Is it worse to believe noise is present when it is absent (e.g. conservative noise model), or ingore sources of noise when they are present?

Ignoring uncertaines that are present:

tbl %>% filter(belief == "HLL", true %in% c("HHL", "HLH", "HHH")) %>% knitr::kable("pandoc")

Conservative: include uncertainty that doesn't exist

tbl %>% filter(true == "HLL", belief %in% c("HHL", "HLH", "HHH")) %>% knitr::kable("pandoc")

We can visualize the whole table (rows are the true scenario, columns are believed scenario). First we just plot ENPV, which justifies our normalization routine, since we see the very strong influence of reality regardless what you believe. Note that we supress grid labels as the units are arbitrary and only relative differences are of interest.

ggplot(tbl) + geom_histogram(aes(ENPV), binwidth=50) + 
  facet_grid(belief ~ true) + 
  theme(axis.text=element_blank())

After normalizing by the optimum that could be achieved for the given scenario, we have the following patterns.

ggplot(tbl) + geom_histogram(aes(normalized_value), binwidth=0.05) + 
  facet_grid(belief ~ true)  + 
  theme(axis.text.y=element_blank()) + 
  scale_x_continuous(breaks=scales::pretty_breaks(n=3)) + 
  coord_cartesian(xlim=c(0, 1.2))

Note that believing that both measurement and implementation error are low (LLL and HLL) gives you consistently the worst outcomes unless both sources of noise truely are low.

ggplot(tbl) + 
  geom_histogram(aes(normalized_value, fill=belief), position="identity", binwidth=.05) + 
  facet_grid(~ true) + 
  theme(axis.text=element_blank())

#write_csv(tbl, "table_files/table.csv")

cboettig/multiple_uncertainty documentation built on May 13, 2019, 2:08 p.m.