R/Myerson.R

In dmi3kno/tidyear: Tidy Economic Analysis (and Decision Analysis) in R

#' Density function for Generalized Lognormal distribution (aka Myerson).
#'
#' @param x vector of quantiles. If `length(n)>1``, the length is taken to be the number required.
#' @param q1 bottom quantile
#' @param q2 middle quantile (50th percentile)
#' @param q3 top quantile
#' @param tl tail percentage
#'
#' @return a vector of probabilities of length equal to `length(x)`.
#' @export
#'
#' @examples
#' dMyerson(0.25, 10, 20, 40)
dMyerson <- function(x,q1,q2,q3,tl=0.5){
  # Probability density function
  br <- (q3-q2)/(q2-q1)

  if(br==1)
    return(dnorm(x, mean=q2, sd=(q3-q2)/qnorm(1-tl/2)))

  q <- qnorm(1-tl/2)*(log((x-q2)*(br-1)/(q3-q2)+1)/log(br))
  num <- qnorm(1-tl/2)*(br-1)
  den <- ((q3-q2)+(x-q2)*(br-1))*log(br)
  num/den*dnorm(q, mean=0, sd=1)
}

#' Cumulative distribution function for Generalized Lognormal distribution (aka Myerson).
#'
#' @param x vector of quantiles. If `length(n)>1``, the length is taken to be the number required.
#' @param q1 bottom quantile
#' @param q2 middle quantile (50th percentile)
#' @param q3 top quantile
#' @param tl tail percentage
#'
#' @return a vector of exceedance probabilities of length equal to `length(x)`.
#' @export
#'
#' @examples
#' pMyerson(0.25, 10, 20, 40)
pMyerson <- function(x,q1,q2,q3,tl=0.5){
  # cumulative density
  br <- (q3-q2)/(q2-q1)

  if(br==1)
    return(pnorm(x, mean=q2, sd=(q3-q2)/qnorm(1-tl/2)))

  nom <- log(1+(x-q2)*(br-1)/(q3-q2))
  den <- log(br)
  q <- qnorm(1-tl/2)*(nom/den)
  pnorm(q, mean=0, sd=1)
}

#' Generate random numbers from Generalized Lognormal distribution (aka Myerson).
#'
#' @param n number of observations. If `length(n)>1``, the length is taken to be the number required.
#' @param q1 bottom quantile
#' @param q2 middle quantile (50th percentile)
#' @param q3 top quantile
#' @param tl tail percentage
#' @param lower lower bound for distribution
#' @param upper upper bound for distribution
#'
#' @return a length `n` vector of random values.
#' @export
#'
#' @examples
#' rMyerson(1, 10, 20, 40)
rMyerson <- function(n,q1,q2,q3,tl=0.5, lower=-Inf, upper=Inf){
  # random number generator
  if(length(n)>1) n <-length(n)
  x <- runif(n)
  norml <- qnorm(x)/qnorm(1-tl/2)
  br <- (q3-q2)/(q2-q1)

  if(br==1)
    return(pmin(pmax(q2+(q3-q2)*norml, lower), upper))

  res <- q2+(q3-q2)*(br^norml-1)/(br-1)
  pmin(pmax(res, lower), upper)
}