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.

> library(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


> 
> #Formula 1
> f <- function(x) log(x)/x
> integrate(f, 1/4, 1)
-0.960906 with absolute error < 3.9e-10
> -log(1/4)^2/2
[1] -0.960906
> 
> #Formula 2 -> D2
> library(gsl)
> f <- function(x) log(1-x)/x
> g <- 3
> b <- 1/4
> integrate(f, b*(g-1)/(g*b-1), (g-1)/(g*b-1))
2.24893 with absolute error < 2.5e-14
> -dilog((g-1)/(g*b-1)) + dilog(b*(g-1)/(g*b-1))
[1] 2.24893
> 
> #Formula 2 -> D1
> g <- 3
> b <- 4
> integrate(f, b*(g-1)/(g*b-1), (g-1)/(g*b-1))
0.7463479 with absolute error < 8.3e-15
> -dilog((g-1)/(g*b-1)) + dilog(b*(g-1)/(g*b-1))
[1] 0.7463479
> 
> #formula 2 -> D2
> 
> f <- function(x) log(1+x*(g-1)/(1-g*b))/x
> g <- 3
> b <- 1/4
> integrate(f, b, 1)
2.24893 with absolute error < 2.5e-14
> -dilog((g-1)/(g*b-1)) + dilog(b*(g-1)/(g*b-1)) 
[1] 2.24893
> 
> #formula 2 -> D2
> 
> f <- function(x) log(1-g*b+x*(g-1))/x
> g <- 3
> b <- 1/4
> integrate(f, b, 1)
0.3271176 with absolute error < 1.2e-10
> -dilog((g-1)/(g*b-1)) + dilog(b*(g-1)/(g*b-1)) - log(1-g*b)*log(b)
[1] 0.3271176
> 
> 
> 
> #Formula 3
> f <- function(x, g, b)
+   log(x)/x/( (g-1)*x+1-g*b )
> 
> I <- function(g, b)
+ {
+   temp <- dilog(b*(g-1)/(g*b-1)) - dilog((g-1)/(g*b-1)) 
+   (log(b)*log(abs(1-b)) - log(b)^2/2 +temp - log(abs(1-g*b))*log(b) )/(1-g*b)
+ }
> 
> b <- 1/4
> integrate(f, b, 1, g=3, b=b)$value
[1] -0.9399057
> I(g=3, b=b)
[1] -0.9399057
> 
> b <- 4
> integrate(f, b, 1, g=3, b=b)$value
[1] 0.1832497
> I(g=3, b=b)
[1] 0.1832497
> 
> 
> #Formula second order moment D3
> f <- function(x)
+     1/(1+(g-1)*sqrt(x))
> 
> I <- function(g)
+     2/(g-1)-2*log(g)/(g-1)^2
> 
> g <- 3.5
> integrate(f, 0, 1)
0.3991158 with absolute error < 8.7e-05
> I(3.5)
[1] 0.3991159
> 
> 
> #formula of regularity conditions - Lemma B1
> 
> f <- function(x)
+   (a+b^x)^m 
> 
> m <- 1; a <- 3; b <- 1/2
> integrate(f, 0, 1)
3.721348 with absolute error < 4.1e-14
> m <- 2; a <- 3; b <- 1/2
> integrate(f, 0, 1)
13.8691 with absolute error < 1.5e-13
> m <- 3; a <- 3; b <- 1/2
> integrate(f, 0, 1)
51.76626 with absolute error < 5.7e-13
> 
> I <- function(m, a, b)
+ {
+   if(m == 0)
+     return(1)
+   if(m > 0)
+   {
+     k <- 1:m
+     return(a^m+ sum(choose(m,k)*a^(m-k)/log(b)*(b^k/k - 1/k)))
+   }else
+   {
+     m <- -m 
+     k <- 1:(m-1)
+     res <- 1/a^{m} - {log({a+b}/{a+1})}/{a^{m}*log(b)}
+     if(m >= 2)
+     res <- res + sum({-1}/{a^k*log(b)*(-m+k)}*({1}/{(a+b)^{m-k}} - {1}/{(a+1)^{m-k}}))
+     return(res)
+   }
+ }
> I(1, a, b)
[1] 3.721348
> I(2, a, b)
[1] 13.8691
> I(3, a, b)
[1] 51.76626
> 
> g <- function(x)
+   (a+b^x)^(-m) 
> m <- 1; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.2691183 with absolute error < 3e-15
> I(-1, a, b)
[1] 0.2691183
> m <- 2; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.07253116 with absolute error < 8.1e-16
> I(-2, a, b)
[1] 0.07253116
> m <- 3; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.01957662 with absolute error < 2.2e-16
> I(-3, a, b)
[1] 0.01957662
> 
> 
> 
> J <- function(m, a, b)
+ {
+   L2ab <- mbbefd:::gendilog(a,b)
+   if(m == 0)
+     return(1)
+   if(m == 1)
+     return(1/(2*a)-(log(a+b) - L2ab)/(a*log(b)))
+   if(m == 2)
+     return(1/a^2/2-log(a+b)/a^2/log(b)+log((a+b)/(a+1))/a^2/log(b)^2 + L2ab/a^2/log(b)-b/log(b)/a^2/(a+b))
+   k <- 1:m
+   J(m-1, a, b)/a - 1/log(b)/a/(-m+1)/(a+b)^{m-1} + I(-m+1, a, b)/log(b)/a/(-m+1)
+ }
> 
> g <- function(x)
+   x/(a+b^x)^(m) 
> 
> m <- 1; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.1375337 with absolute error < 1.5e-15
> J(1, a, b)
[1] 0.1375337
> m <- 2; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.0378636 with absolute error < 4.2e-16
> J(2, a, b)
[1] 0.0378636
> m <- 3; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.01043275 with absolute error < 1.2e-16
> J(3, a, b)
[1] 0.01043275
> m <- 4; a <- 3; b <- 1/2
> integrate(g, 0, 1)
0.002876944 with absolute error < 3.2e-17
> J(4, a, b)
[1] 0.002876944
> 
> 
> #formula of regularity conditions - Lemma B2
> f <- function(x)
+   b^x*log(b)/((a+b^x)^m)
> 
> I <- function(m, a, b)
+ {
+   if(m == 1)
+     return(log((a+b)/(a+1)))
+   1/(-m+1)*(1/(a+b)^{m-1} - 1/(a+1)^{m-1})
+ }
> 
> m <- 1; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.1335314 with absolute error < 1.5e-15
> I(1, a, b)
[1] -0.1335314
> 
> m <- 2; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.03571429 with absolute error < 4e-16
> I(2, a, b)
[1] -0.03571429
> 
> 
> #formula of regularity conditions - Lemma B3
> 
> 
> f <- function(x)
+   x*b^x*log(b)/(a+b^x)^m
> 
> 
> 
> I <- function(m, a, b)
+ {
+   if(m == 1)
+     return(log(a+b)-mbbefd:::gendilog(a,b))
+   res <- b/a/(a+b) - log((a+b)/(a+1))/a/log(b)
+   if(m>2)
+   {
+     l <- 1:(m-2)
+     res <- res - sum(1/a^l/log(b)/(m-1)/(-m+1+l)*(1/(a+b)^(m-1-l) - 1/(a+1)^(m-1-l)) ) 
+   }
+   res  
+ }
>   
> m <- 1; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.06058021 with absolute error < 6.7e-16
> I(m, a,b)
[1] -0.06058021
> 
> m <- 2; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.01659598 with absolute error < 1.8e-16
> I(m, a,b)
[1] -0.01659598
> 
> m <- 3; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.004550746 with absolute error < 5.1e-17
> I(m, a,b)
[1] -0.02518345
> 
> m <- 4; a <- 3; b <- 1/2
> integrate(f, 0, 1)
-0.001248997 with absolute error < 1.4e-17
> I(m, a,b)
[1] -0.02003778
> 
> 
> #formula of regularity conditions - Lemma B4
> 
> f <- function(x)
+   x^2*b^x*log(b)/(a+b^x)^m
> 
> I <- function(m, a, b)
+ {
+   L2ab <- mbbefd:::gendilog(a,b)
+   
+   if(m == 2)
+     res <- -1/(a+b)+1/a+2/a/log(b)*(-log(a+b)+L2ab)
+   else if(m == 3)
+     res <- -1/2/(a+b)^2+1/a^2/2-log(a+b)/a^2/log(b)+log((a+b)/(a+1))/a^2/log(b)^2 + L2ab/a^2/log(b)-b/log(b)/a^2/(a+b)
+   else
+     res <- NA
+   res  
+ }
> 
> m <- 2; a <- 3.1; b <- 1/2
> integrate(f, 0, 1)
-0.01007996 with absolute error < 1.1e-16
> I(m, a,b)
[1] -0.01007996
> 
> m <- 3; a <- 3.1; b <- 1/2
> integrate(f, 0, 1)
-0.002719957 with absolute error < 3e-17
> I(m, a,b)
[1] -0.002719957
> 
> proc.time()
   user  system elapsed 
   1.26    0.12    1.37