R Under development (unstable) (2023-08-27 r85021 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #test integral for second order moment
> 
> f <- function(x, a, b)
+ {
+   if(a == 0)
+     return(x*log(b))
+   if(b == 1)
+     return(rep(log(a+1), length(x)))
+   log(a+b^x)
+ }
> 
> 
> a <- -1/2
> sapply(1:3+1/2, function(b)
+   c(a=a, b=b, check=integrate(f, 0, 1, a=a, b=b)$value, th=mbbefd:::gendilog(a, b))
+ )
            [,1]        [,2]       [,3]
a     -0.5000000 -0.50000000 -0.5000000
b      1.5000000  2.50000000  3.5000000
check -0.3300563  0.05299366  0.2800041
th    -0.3300563  0.05299366  0.2800041
> 
> a <- 0
> sapply(1:3, function(b)
+   c(a=a, b=b, check=integrate(f, 0, 1, a=a, b=b)$value, th=mbbefd:::gendilog(a, b))
+ )
      [,1]      [,2]      [,3]
a        0 0.0000000 0.0000000
b        1 2.0000000 3.0000000
check    0 0.3465736 0.5493061
th       0 0.3465736 0.5493061
> 
> a <- 1
>   sapply(1:3, function(b)
+ c(a=a, b=b, check=integrate(f, 0, 1, a=a, b=b)$value, th=mbbefd:::gendilog(a, b, checkparam=FALSE))
+ )
           [,1]      [,2]     [,3]
a     1.0000000 1.0000000 1.000000
b     1.0000000 2.0000000 3.000000
check 0.6931472 0.8862177 1.016654
th    0.6931472 0.8862177 1.016654
> 
> a <- 2
> sapply(1:3, function(b)
+   c(a=a, b=b, check=integrate(f, 0, 1, a=a, b=b)$value, th=mbbefd:::gendilog(a, b, checkparam=FALSE))
+ )
          [,1]     [,2]     [,3]
a     2.000000 2.000000 2.000000
b     1.000000 2.000000 3.000000
check 1.098612 1.232791 1.329374
th    1.098612 1.232791 1.329374
> 
> a <- 3
> sapply(1:3, function(b)
+   c(a=a, b=b, check=integrate(f, 0, 1, a=a, b=b)$value, th=mbbefd:::gendilog(a, b, checkparam=FALSE))
+ )
          [,1]     [,2]    [,3]
a     3.000000 3.000000 3.00000
b     1.000000 2.000000 3.00000
check 1.386294 1.489181 1.56596
th    1.386294 1.489181 1.56596
> 
> 
> vartheo <- function(a, b)
+ {
+   EX <- mmbbefd(1, a, b)
+   
+   mmbbefd(2, a, b) - EX^2
+ }
> 
> require(mbbefd)
Loading required package: mbbefd
Loading required package: fitdistrplus
Loading required package: MASS
Loading required package: survival
Loading required package: alabama
Loading required package: numDeriv
Loading required package: Rcpp
Package:  mbbefd
Version:  0.8.11
Date:     
BugReport: https://github.com/spedygiorgio/mbbefd/issues


> n <- 1e4
> x <- rmbbefd(n, 2, 1/2)
> c(var(x), vartheo(2, 1/2))
[1] 0.09946506 0.09966670
> 
> x <- rmbbefd(n, -1/2, 3)
> c(var(x), vartheo(-1/2, 3))
[1] 0.1366748 0.1355999
> 
> x <- rmbbefd(n, -1/2, 2)
> c(var(x), vartheo(-1/2, 2))
[1] 0.1105597 0.1117094
> 
> x <- rmbbefd(n, 2, 1/5)
> c(var(x), vartheo(2, 1/5))
[1] 0.1222539 0.1230419
> 
> 
> proc.time()
   user  system elapsed 
   1.04    0.12    1.15