[FEATURE] Add composition scaling class for Bernstein Functions

[hsloot]

Summary

Compositions of Bernstein functions are again Bernstein functions. This holds in particular for the inner functions being a linear function. This is easy to implement and will allow a larger variety of Bernstein functions to be created.

Proposal

A possible implementation

## allClass-S4.R

CompositScaledBernsteinFunction <- setClass("CompositScaledBernsteinFunction",
  contains  = "BernsteinFunction",
  slots = c(cscale = "numeric", original = "BernsteinFunction"))

## allValidity-S4.R
setValidity("CompositScaledBernsteinFunction",
  function(object) {
    qassert(object@cscale, "N1(0,)")

    invisible(TRUE)
  })

## allInitialize-S4.R

setMethod("initialize", "CompositScaledBernsteinFunction",
  function(.Object, cscale = 1, original = AlgebraicBernsteinFunction2()) { # nolint
    .Object@cscale <- cscale
    .Object@original <- original
    validObject(.Object)

    invisible(.Object)
  })

## valueOf-S4.R

setMethod("valueOf", "CompositScaledBernsteinFunction",
  function(object, x, difference_order = 0L, cscale = 1, n = 1, k = 0, ...,
           method = c("default", "levy", "stieltjes"),
           tolerance = .Machine$double.eps^0.5) {
    method <- match.arg(method)
    if (!"default" == method) {
      return(callNextMethod())
    }
    assert(combine = "or",
           check_numeric(x, lower = 0, min.len = 1L, any.missing = FALSE),
           check_complex(x, min.len = 1L, any.missing = FALSE))
    qassert(difference_order, "X1[0,)")
    qassert(cscale, "N1(0,)")
    qassert(n, "X1(0,)")
    qassert(k, "N1[0,)")

    valueOf(object@original, x, difference_order, cscale*object@cscale, n, k, ...)
  })

As a showcase, the following shows how this can be used to create the InverseGaussianBernsteinFunction based on the Algebraic Bernstein function no. 2.

## allClass-S4.R

AlgebraicBernsteinFunction2 <- setClass("AlgebraicBernsteinFunction2",
  contains = "CompleteBernsteinFunction",
  slots = c(alpha = "numeric"))

## allValidity-S4.R

setValidity("AlgebraicBernsteinFunction2",
  function(object) {
    qassert(object@alpha, "N1(0,)")

    invisible(TRUE)
  })

## allInitialize-S4.R

setMethod("initialize", "AlgebraicBernsteinFunction2",
  function(.Object, alpha = log2(2 - 0.5)) { # nolint
    .Object@alpha <- alpha
    validObject(.Object)

    invisible(.Object)
  })

## levyDensity-S4.R

setMethod("levyDensity", "AlgebraicBernsteinFunction2",
  function(object) {
    structure(
      function(x) {
        object@alpha / gamma(1 - object@alpha) * x ^ (-1 - object@alpha) * exp(-x)
      },
      lower = 0, upper = Inf, type = "continuous"
    )
  })

## stieltjesDensity-S4.R

setMethod("stieltjesDensity", "AlgebraicBernsteinFunction2",
  function(object) {
    structure(
      function(x) {
        sin(object@alpha * pi) / pi * (x - 1) ^ (object@alpha - 1) / x
      },
      lower = 1, upper = Inf, type = "continuous"
    )
  })

eta <- 2
valueOf(InverseGaussianBernsteinFunction(eta = eta), seq(0, 2, by = 0.25), 0L, cscale = 1)

valueOf(
  ScaledBernsteinFunction(
    scale = eta,
    original = CompositScaledBernsteinFunction(
      cscale = 2 / eta^2,
      original = AlgebraicBernsteinFunction2(
        alpha = 0.5
        )
      )
    ),
  seq(0, 2, by = 0.25))