R/cf_LogRV_FisherSnedecorNC.R

In gajdosandrej/CharFunToolR: Numerical Computation Cumulative Distribution Function and Probability Density Function from Characteristic Function

#' @title Characteristic function of a linear combinationof independent
#' LOG-TRANSFORMED non-central Fisher-Snedecor random variables
#'
#' @description
#' \code{cf_LogRV_FisherSnedecorNC(t, df1, df2, delta, coef, niid, tol)} evaluates characteristic function
#' of a linear combination (resp. convolution) of independent LOG-TRANSFORMED non-central Fisher-Snedecor
#' random variables, with distributions \eqn{F(df1_i,df2_i,\delta_i)}.
#'
#' That is, \code{cf_LogRV_FisherSnedecorNC} evaluates the characteristic function
#' \eqn{cf(t)} of  \eqn{Y = coef_i*log(X_1) +...+ coef_N*log(X_N)}, where \eqn{X_i ~ F(df1_i,df2_i,\delta_i)}
#' are inedependent RVs, with \eqn{df1_i} and \eqn{df2_i} degrees of freedom, and the noncentrality parameters \eqn{delta_i >0},
#' for \eqn{i =  1,...,N}.
#'
#' The characteristic function of \eqn{Y = log(X) with X ~ F(df1,df2,\delta)} is Poisson mixture
#' of the CFs of the shifted log-transformed central F RVs of the form
#' \deqn{cf(t) = cf_LogRV_FisherSnedecorNC(t,df1,df2,\delta) = exp(-\delta/2) sum_{j=1}^Inf (\delta/2)^j/j! *
#' exp(1i*t*(df1+2*j)/df1) * cf_LogRV_FisherSnedecor(t,df1+2*j,df2),}
#' where cf_LogRV_FisherSnedecor(t,df1,df2) denotes CF of log-transformed centrally distributed
#' F RVs with parameters df1 and df2. For more details on  the non-central
#' F distribution see \code{cf_FisherSnedecorNC}.
#' Alternatively,
#' \deqn{cf(t) = (df2/df1)^(1i*t) * gamma(df1/2 + 1i*t) / gamma(df1/2) * gamma(df2/2 - 1i*t) / gamma(df2/2) * 1F1(-1i*t;df1/2;-delta/2),}
#' where \eqn{1F1(a;b;z)} is the confluent hypergeometric function, also known as the Kummer function \eqn{M(a,b,z)}.
#' Hence,the characteristic function of \eqn{Y  = coef(1)*Y1 + ... + coef(N)*YN}
#' is  \eqn{cf_Y(t) =  cf_Y1(coef(1)*t) * ... * cf_YN(coef(N)*t)}, where \eqn{cf_Yi(t)}
#' is evaluated with the parameters \eqn{df1_i}, \eqn{df2_i}, and \eqn{delta_i}.
#'
#' @param t vector or array of real values, where the CF is evaluated.
#' @param df1 vector of the  degrees of freedom \code{df1 > 0}. If empty, default value is \code{df1 = 1}.
#' @param df2 vector of the  degrees of freedom \code{df2 > 0}. If empty, default value is \code{df2 = 1}.
#' @param delta vector of non-centrality parameters.
#' @param coef vector of the coefficients of the linear combination of the Beta distributed random variables.
#' If coef is scalar, it is assumed that all coefficients are equal. If empty, default value is \code{coef = 1}.
#' @param niid scalar convolution coeficient \code{niid}, such that \eqn{Z = Y + ... + Y}
#' is sum of \eqn{niid} iid random variables \eqn{Y}, where each \eqn{Y = sum_{i=1}^N coef(i) * log(X_i)}
#' is independently and identically distributed random variable. If empty, default value is \code{niid = 1}.
#' @param tol tolerance factor for selecting the Poisson weights, i.e. such that \eqn{PoissProb > tol}.
#' If empty, default value is \code{tol = 1e-12}.
#'
#' @return Characteristic function \eqn{cf(t)} of a linear combination
#' of independent LOG-TRANSFORMED non-central Fisher-Snedecor random variables.
#'
#' @seealso For more details see WIKIPEDIA:
#' \url{https://en.wikipedia.org/wiki/Noncentral_F-distribution}.
#'
#' @family Continuous Probability Distribution
#' @family Non-central Probability Distribution
#'
#' @note Ver.: 20-Sep-2018 19:44:50 (consistent with Matlab CharFunTool v1.3.0, 10-Aug-2018 15:46:49).
#'
#' @example R/Examples/example_cf_LogRV_FisherSnedecorNC.R
#'
#' @export
#'
cf_LogRV_FisherSnedecorNC <- function(t, df1, df2, delta, coef, niid, tol) {
  ## CHECK THE INPUT PARAMETERS
  if(missing(df1)) {
    df1 <- vector()
  }
  if(missing(df2)) {
    df2 <- vector()
  }
  if(missing(delta)) {
    delta <- vector()
  }
  if(missing(coef)) {
    coef <- vector()
  }
  if(missing(niid)) {
    niid <- numeric()
  }
  if(missing(tol)) {
    tol <- numeric()
  }

  if(length(df1) == 0) {
    df1 <- 1
  }
  if(length(df2) == 0) {
    df2 <- 1
  }
  if(length(delta) == 0) {
    delta <- 0
  }
  if(length(coef) == 0) {
    coef <- 1
  }
  if(length(niid) == 0) {
    niid <- 1
  }
  if(length(tol) == 0) {
    tol <- 1e-12
  }

  ## SET THE COMMON SIZE of the parameters
  l_max <- max(c(length(df1), length(df2), length(delta), length(coef)))
  if (l_max > 1) {
    if (length(df1) == 1) {
      df1 <- rep(df1, l_max)
    }
    if (length(df2) == 1) {
      df2 <- rep(df2, l_max)
    }
    if (length(delta) == 1) {
      delta <- rep(delta, l_max)
    }
    if (length(coef) == 1) {
      coef <- rep(coef, l_max)
    }
    if ((any(lengths(list(
      coef, df1, df2, delta
    )) < l_max))) {
      stop("Input size mismatch.")
    }
  }

  ## Characteristic function of a linear combination of independent nc F RVs
  szt <- dim(t)
  t <- c(t)
  cf <- rep(1, length(t))
  for(i in 1:length(coef)) {
    cf <- cf * cf_ncLogRVFisherSnedecor(coef[i] * t, df1[i], df2[i], delta[i], tol)
  }
  dim(cf) <- szt
  cf[t==0] <- 1

  if(length(niid) > 0) {
    if(length(niid) == 1) {
      cf <- cf ^ niid
    } else {
      stop("niid should be a scalar (positive integer) value")
    }
  }

  return(cf)
}

## Function funCF
cf_ncLogRVFisherSnedecor <- function(t, df1, df2, delta, tol) {
  # cf_ncLogRVFisherSnedecor Characteristic function of the distribution of
  # the of the distribution of the log-transformed non-central
  # Fisher-Snedecor RV with df1 and df2 degrees of freedom and the
  # non-centrality parameter delta > 0.

  f <- 0
  delta  <- delta/2
  if(delta == 0) {
    # Deal with the central distribution
    f <- cf_LogRV_FisherSnedecor(t, df1, df2)
  } else if(delta > 0) {
    # Sum the Poisson series of CFs of independent log-transformed F RVs,
    # poisspdf(j,delta).*cf_LogRV_FisherSnedecor(t*(df1+2*j)/df1,df1+2*j,df2)
    j0 <- floor(delta / 2)
    p0 <- exp(-delta + j0 * log(delta) - log(gamma(j0 + 1)))
    f <- f + p0 * exp(1i * t * log((df1 + 2 * j0) / df1)) * cf_LogRV_FisherSnedecor(t, df1 + 2 * j0, df2)
    p <- p0
    j <- j0 - 1
    while(j >= 0 && p > tol) {
      p <- p * (j + 1) / delta
      f <- f + p * exp(1i*t*log((df1+2*j)/df1)) * cf_LogRV_FisherSnedecor(t, df1 + 2 * j, df2)
      j <- j - 1
    }
    p <- p0
    j <- j0 + 1
    i <- 0
    while(p > tol && i <= 5000) {
      p <- p * delta / j
      f <- f + p * exp(1i*t*log((df1+2*j)/df1)) * cf_LogRV_FisherSnedecor(t, df1 + 2 * j, df2)
      j <- j + 1
      i <- i + 1
    }
    if(i == 5000) {
      warning("No convergence.")
    }
    return (f)
  } else {
    stop("delta should be nonnegative.")
  }
}