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


> 
> #test of MBBEFD(a,b) distribution
> n <- 1e4
> 
> set.seed(567)
> x <- rmbbefd(n, 2, 1/2)
> y <- rmbbefd(n, -1/2, 2)
> xD4 <- rmbbefd(n, Inf, 1/3)
> xD3 <- rmbbefd(n, 0, 1/3)
> xD5 <- rmbbefd(n, -1, 3)
> 
> #test CDF
> z <- 0:8/8
> cbind(ecdf(x)(z), pmbbefd(z, 2, 1/2))
        [,1]       [,2]
 [1,] 0.0000 0.00000000
 [2,] 0.0578 0.05690493
 [3,] 0.1149 0.11200942
 [4,] 0.1713 0.16520092
 [5,] 0.2198 0.21638838
 [6,] 0.2660 0.26550189
 [7,] 0.3100 0.31249201
 [8,] 0.3555 0.35732870
 [9,] 1.0000 1.00000000
> 
> cbind(ecdf(y)(z), pmbbefd(z, -1/2, 2))
        [,1]       [,2]
 [1,] 0.0000 0.00000000
 [2,] 0.0779 0.07663552
 [3,] 0.1386 0.13726434
 [4,] 0.1833 0.18626055
 [5,] 0.2227 0.22654092
 [6,] 0.2549 0.26012531
 [7,] 0.2841 0.28845700
 [8,] 0.3118 0.31259484
 [9,] 1.0000 1.00000000
> 
> #test EC
> cbind(eecf(x)(z), ecmbbefd(z, 2, 1/2))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.1538975 0.1538776
 [3,] 0.2986830 0.2988821
 [4,] 0.4345831 0.4353077
 [5,] 0.5622590 0.5634651
 [6,] 0.6822638 0.6836774
 [7,] 0.7951934 0.7962775
 [8,] 0.9010829 0.9016043
 [9,] 1.0000000 1.0000000
> 
> cbind(eecf(y)(z), ecmbbefd(z, -1/2, 2))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.1508709 0.1514407
 [3,] 0.2910350 0.2921265
 [4,] 0.4230304 0.4242127
 [5,] 0.5484774 0.5492894
 [6,] 0.6682255 0.6685629
 [7,] 0.7829882 0.7829694
 [8,] 0.8934175 0.8932498
 [9,] 1.0000000 1.0000000
> 
> #test mean
> mean(x)
[1] 0.7878581
> mmbbefd(1, 2, 1/2)
[1] 0.7891032
> 
> mean(y)
[1] 0.7940801
> mmbbefd(1, -1/2, 2)
[1] 0.7924813
> 
> mean(xD4)
[1] 0.6078659
> (1/3-1)/log(1/3)
[1] 0.6068262
> 
> mean(xD3)
[1] 1
> mean(xD5)
[1] NaN
> 
> #second order moment
> mean(x^2)
[1] 0.722214
> mmbbefd(2, 2, 1/2)
[1] 0.7223506
> 
> mean(y^2)
[1] 0.7419741
> mmbbefd(2, -1/2, 2)
[1] 0.739736
> 
> mean(xD4^2)
[1] 0.4981949
> mmbbefd(2, Inf, 1/3)
[1] 0.4978878
> 
> mean(xD3^2)
[1] 1
> mmbbefd(2, 0, 1/3)
[1] 1
> 
> mean(xD5^2)
[1] NaN
> 
> #total loss
> etl(x)
[1] 0.6023
> tlmbbefd(2, 1/2)
[1] 0.6
> 
> 
> etl(y)
[1] 0.6678
> tlmbbefd(-1/2, 2)
[1] 0.6666667
> 
> 
> #test quantile
> cbind(quantile(y, probs=0:10/10), qmbbefd(0:10/10, -1/2, 2))
            [,1]      [,2]
0%   0.000608317 0.0000000
10%  0.168459642 0.1699250
20%  0.430643652 0.4150375
30%  0.819179690 0.8073549
40%  1.000000000 1.0000000
50%  1.000000000 1.0000000
60%  1.000000000 1.0000000
70%  1.000000000 1.0000000
80%  1.000000000 1.0000000
90%  1.000000000 1.0000000
100% 1.000000000 1.0000000
> qmbbefd(1/2, -1/2, 2)
[1] 1
> 
> 
> z <- seq(0, 1, length=101)
> plot(z, pmbbefd(z, -1/2, 2), type="l", ylim=c(0, 1-tlmbbefd(-1/2, 2)))
> 
> plot(z, qmbbefd(z, -1/2, 2), type="l", xlim=c(0, 1-tlmbbefd(-1/2, 2)))
> 
> 
> 
> #test density
> 
> z <- sort(c(1, seq(-0.1,1.1, length=101)))
> plot(density(x), ylim=c(0,1))
> lines(z, dmbbefd(z, 2, 1/2), col="red")
> 
> 
> plot(density(y), ylim=c(0,1))
> lines(z, dmbbefd(z, -1/2, 2), col="red")
> 
> proc.time()
   user  system elapsed 
   1.32    0.17    1.50