R/distribution-seasonal-temperatures.R

In EarthSystemDiagnostics/ecustools: A miscellaneous collection of R functions for working with time-series and paleo data.

#' Distribution of temperatures in seasonal cycle
#'
#' @description Probability density function and random generation for seasonal
#'   temperatures
#' @details The seasonal temperature cycle is approximated as a sine wave with
#'   gaussian noise. The probability density function is obtained by numerically
#'   convoluting a scaled sine wave and normal PDF with mean = 0. Random
#'   deviates are generated by adding gaussian noise to points sampled uniformly
#'   from the scaled sine wave on the interval 0 to 2*pi
#' @param x vector of temperatures
#' @param mean.t mean annual temperature
#' @param amp.t amplitude of the seasonal cycle
#' @param sd.t inter-annual standard deviation of temperature at a given point
#'   in the seasonal cycle
#' @param return return densities for temperatures x, or a function to do so,
#'   Default: c("density", "FUN")
#' @param res resolution of the temperature vectors to be convolved, Default:
#'   0.01
#' @name SeasonalCycle
#' @author Andrew Dolman
#' @examples
#' \dontrun{
#' if(interactive()){
#' amp.t <- 10
#' mean.t <- 25
#' sd.t <- 0.75
#'
#' x <- seq(mean.t-amp.t, mean.t+amp.t, length.out = 1000)
#' Z <- rSeas(n = 100000, mean.t, amp.t, sd.t)
#'
#' hist(Z, freq=F,breaks=100)
#' lines(x, dSeas(x, mean.t, amp.t, sd.t, res = 0.1), col = "Blue")
#' lines(x, dSeas(x, mean.t, amp.t, sd.t, res = 0.01), col = "Green")
#' lines(x, dSeas(x, mean.t, amp.t, sd.t, res = 0.001), col = "Red")
#'  }
#'}
NULL

#' @rdname SeasonalCycle
#' @export
dSeas <- function(x=NULL, mean.t, amp.t, sd.t, return = c("density", "FUN"), res = 0.01) {
  function_deprecated("stattools")
  rng <- max(c(amp.t, sd.t * 5))
  lnth <- (2 * rng) / res

  # make vector length with many factors for fast FFT
  lnth2 <- (2 ^ ceiling(log(lnth, base = 2)))

  if (lnth2 > 5000)  warning("Sequences to be convolved are very long,
                             decrease the resolution or rescale mean.t and amp.t")

  # calculate new resolution for scaling
  res2 <- (2 * rng) / lnth2

  xin <- seq(mean.t - rng, mean.t + rng, length.out = lnth2)
  zin <- dSin(xin, mean.t, amp.t)
  #plot(xin, zin, type = "l")

  yin <- dnorm(xin, mean.t, sd.t)
  #plot(xin, yin, type = "l")

  cv <- convolve(zin, yin, type = "open")
  cv <- cv / (sum(cv) * res2)

  x2 <- seq(mean.t - 2 * rng, mean.t + 2 * rng, length.out = length(cv))

  FUN <- approxfun(x2, cv, yleft = 0, yright = 0)

  ret <- match.arg(return)

  switch(ret,
         density = return(FUN(x)),
         FUN = return(FUN))
}


#' @rdname SeasonalCycle
#' @param n number of samples
#' @export
rSeas <- function(n, mean.t, amp.t, sd.t){
  function_deprecated("stattools")
  x <- runif(n, 0, 2*pi)
  y <- mean.t + sin(x) * amp.t/2
  z <- y + rnorm(n, 0, sd.t)
  return(z)
}



#' Density of a scaled sin wave
#' @param mean.t mean of sine wave
#' @param amp.t amplitude of sine wave
#' @export
dSin <- function(x, mean.t, amp.t) {
  function_deprecated("stattools")
  d <- rep(0, length(x))

  # Do not evaluate where known to be NaN
  ind <- x > (mean.t - amp.t / 2) & x < (mean.t + amp.t / 2)
  d[ind] <- 1 / (pi * sqrt(1 - ((2 * x[ind]) / amp.t - (2 * mean.t) / amp.t) ^ 2))

  2 * d / amp.t
}