[FEAT] Improve Bernstein Function representation

[hsloot]

Summary

The current (public) representations of Bernstein functions are not unique; thus potentially leading to the issue that identical(bf1, bf2) is false even if bf1and bf2 represent the same Bernstein function. Additionally, the current system leads to very confusing printed representations for larger sums; see reprex below.

Proposal

The public API should have a one-for-all Bernstein function representation that allows representing sums, scaling, and composite scaling. The representation should be unique, i.e., a Bernstein function must not have more than one representation (as far as this is possible).

Additional context

library(rmo)

bf1 <- AlphaStableBernsteinFunction(alpha = 0.4)
bf2 <- PoissonBernsteinFunction(eta = 2)
bf3 <- ExponentialBernsteinFunction(lambda = 0.5)
bf4 <- ParetoBernsteinFunction(alpha = 0.2, x0 = 1e-2)
bf5 <- GammaBernsteinFunction(a = 2)
bf6 <- InverseGaussianBernsteinFunction(eta = 0.3)

bf7 <- SumOfBernsteinFunctions(bf1, bf2)
bf8 <- SumOfBernsteinFunctions(bf2, bf1)
bf7
#> An object of class "SumOfBernsteinFunctions"
#> - first:
#>  An object of class "AlphaStableBernsteinFunction"
#>  - alpha: 0.4
#> - second:
#>  An object of class "PoissonBernsteinFunction"
#>  - eta: 2
bf8
#> An object of class "SumOfBernsteinFunctions"
#> - first:
#>  An object of class "PoissonBernsteinFunction"
#>  - eta: 2
#> - second:
#>  An object of class "AlphaStableBernsteinFunction"
#>  - alpha: 0.4
identical(bf7, bf8)
#> [1] FALSE

bf9 <- SumOfBernsteinFunctions(
    ScaledBernsteinFunction(0.5, bf1),
    ScaledBernsteinFunction(0.5, bf2)
)
bf10 <- ScaledBernsteinFunction(
    0.5,
    SumOfBernsteinFunctions(bf1, bf2)
)
bf9
#> An object of class "SumOfBernsteinFunctions"
#> - first:
#>  An object of class "ScaledBernsteinFunction"
#>  - scale: 0.5
#>  - original:
#>      An object of class "AlphaStableBernsteinFunction"
#>      - alpha: 0.4
#> - second:
#>  An object of class "ScaledBernsteinFunction"
#>  - scale: 0.5
#>  - original:
#>      An object of class "PoissonBernsteinFunction"
#>      - eta: 2
bf10
#> An object of class "ScaledBernsteinFunction"
#> - scale: 0.5
#> - original:
#>  An object of class "SumOfBernsteinFunctions"
#>  - first:
#>      An object of class "AlphaStableBernsteinFunction"
#>      - alpha: 0.4
#>  - second:
#>      An object of class "PoissonBernsteinFunction"
#>      - eta: 2
identical(bf7, bf8)
#> [1] FALSE

bf11 <- SumOfBernsteinFunctions(
    bf1,
    SumOfBernsteinFunctions(
        bf2,
        SumOfBernsteinFunctions(
            bf3,
            SumOfBernsteinFunctions(
                bf4,
                SumOfBernsteinFunctions(
                    bf5,
                    bf6
                )
            )
        )
    )
)
bf11
#> An object of class "SumOfBernsteinFunctions"
#> - first:
#>  An object of class "AlphaStableBernsteinFunction"
#>  - alpha: 0.4
#> - second:
#>  An object of class "SumOfBernsteinFunctions"
#>  - first:
#>      An object of class "PoissonBernsteinFunction"
#>      - eta: 2
#>  - second:
#>      An object of class "SumOfBernsteinFunctions"
#>      - first:
#>          An object of class "ExponentialBernsteinFunction"
#>          - lambda: 0.5
#>      - second:
#>          An object of class "SumOfBernsteinFunctions"
#>          - first:
#>              An object of class "ParetoBernsteinFunction"
#>              - alpha: 0.2
#>              - x0: 0.01
#>          - second:
#>              An object of class "SumOfBernsteinFunctions"
#>              - first:
#>                  An object of class "GammaBernsteinFunction"
#>                  - a: 2
#>              - second:
#>                  An object of class "InverseGaussianBernsteinFunction"
#>                  - eta: 0.3

^{Created on 2023-02-04 with reprex v2.0.2}

[hsloot]

Consider representing Bernstein functions as a vector space over a basis consisting of one to two parametric Bernstein functions. For example

setClass(
    "VectorBernsteinFunction",
    contains = "BernsteinFunction",
    slots = c(
        coefficients = "numeric",
        basis = "list"              ## list of basis Bernstein functions
    )

setValidity(
    "VectorBernsteinFunction",
    function(object) {
        if (!qtest(object@coefficients, "N+(0,)")) {
            return(error_msg_domain("coefficients", "N+(0,)"))
        }
        if (!test_list(object@basis, types = "BasisBernsteinFunction", any.missing = FALSE, len = length(coefficients))) {
            return(error_msg_domain("basis")
        }

        invisible(TRUE)
    })