In EarthSystemDiagnostics/spectraluncertainty:

knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE, cache = TRUE,
  fig.width = 7, fig.height = 4.5, fig.pos = "H", dpi = 300, autodep = TRUE,
  dev = "cairo_pdf"
)

# detach("package:spectraluncertainty", unload = TRUE)
# detach("package:sedproxy", unload = TRUE)
library(spectraluncertainty)
library(sedproxy)
library(tidyverse)

GetCoverage <- function(spec.pars, tau_smooth = spec.pars$tau_r){


  stopifnot(tau_smooth %% spec.pars$delta_t == 0)

  sed.pars <- SedproxyParsFromSpecPars(spec.pars)
  pfm <- do.call(ClimToProxyClim, sed.pars)

  # smooth pfm
  flt.wdth <- tau_smooth / spec.pars$delta_t 
  flt <- rep(1/flt.wdth, flt.wdth)

  pfm$simulated.proxy <- as.data.frame(
    apply(pfm$simulated.proxy, 2, function(x) stats::filter(x, flt))
    )

pfm.err <- GetErrorComponents(pfm) %>% 
  select(timepoints, Bioturbation, Measurement, Aliasing.M, Aliasing.Y, Aliasing.YM) %>% 
  rename(Indep.error = Measurement,
         Aliasing.seasonal = Aliasing.M,
         Aliasing.stochastic = Aliasing.Y,
         Aliasing.total = Aliasing.YM) %>% 
  gather(component, error, -timepoints)

pes <- do.call(ProxyErrorSpectrum, spec.pars)
pv <- IntegrateErrorSpectra(pes, tau_smooth = tau_smooth)
pe <- GetProxyError(pv, timescale = tau_smooth) 

sum.as.var <- function(x){sqrt(sum(x^2))}

# combine the two independent error components
ind.err <- pe %>% 
  filter(component %in% c("Meas.error", "Individual.variation")) %>% 
  mutate(component = "Indep.error") %>% 
  group_by(component, smoothed.resolution) %>% 
  summarise_if(is.numeric, sum.as.var) 

# combine the two aliasing components
alias.tot <- pe %>% 
  filter(component %in% c("Aliasing.stochastic", "Aliasing.seasonal")) %>% 
  mutate(component = "Aliasing.total") %>% 
  group_by(component, smoothed.resolution) %>% 
  summarise_if(is.numeric, sum.as.var)


pe <- bind_rows(pe, ind.err, alias.tot)

covrg <- left_join(pfm.err, pe) %>% 
  group_by(component) %>% 
  mutate(covered = abs(error) <= inc.f.zero) %>% 
  filter(complete.cases(timepoints)) %>% 
  summarise(coverage = sum(covered) / n()) 

return(covrg)
}

spec.pars <- GetSpecPars("Mg_Ca", T = 1e04 + 100, sig.sq_a = 0, tau_p = 0.5,
                         sigma.meas = 0.1, sigma.ind = 2, N = 30)


coverage <- tibble(rep = 1:30, tau_smooth = 1000) %>% 
  group_by(rep) %>% 
  do({
    suppressMessages(GetCoverage(spec.pars, tau_smooth = .$tau_smooth))
  })

coverage %>% 
  ggplot(aes(x = component, y = coverage)) +
  geom_violin() +
  geom_point(position = position_jitter(width = 0.2)) +
  geom_hline(yintercept = 0.6827) +
  expand_limits(y = c(0, 1))

EarthSystemDiagnostics/spectraluncertainty documentation built on Jan. 8, 2020, 9:43 a.m.