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/5)
> y <- rMBBEFD(n, 3, 2)
> x3 <- rMBBEFD(n, 3, 1)
> x4 <- rMBBEFD(n, 3, 1/3)
> x5 <- rMBBEFD(n, 3, 0)
> 
> #test CDF
> z <- 0:8/8
> cbind(ecdf(x)(z), round(pMBBEFD(z, 2, 1/5), 6))
        [,1]     [,2]
 [1,] 0.0000 0.000000
 [2,] 0.0543 0.052771
 [3,] 0.1134 0.110191
 [4,] 0.1761 0.171599
 [5,] 0.2393 0.236068
 [6,] 0.3009 0.302451
 [7,] 0.3674 0.369453
 [8,] 0.4328 0.435732
 [9,] 1.0000 1.000000
> 
> cbind(ecdf(y)(z), round(pMBBEFD(z, 3, 2), 6))
        [,1]     [,2]
 [1,] 0.0000 0.000000
 [2,] 0.2456 0.249240
 [3,] 0.3895 0.388908
 [4,] 0.4758 0.477964
 [5,] 0.5392 0.539504
 [6,] 0.5869 0.584428
 [7,] 0.6216 0.618551
 [8,] 0.6459 0.645262
 [9,] 1.0000 1.000000
> 
> #test EC
> cbind(eecf(x)(z), round(ecMBBEFD(z, 2, 1/5), 6))
           [,1]     [,2]
 [1,] 0.0000000 0.000000
 [2,] 0.1603038 0.160391
 [3,] 0.3113117 0.311703
 [4,] 0.4521929 0.453220
 [5,] 0.5827969 0.584359
 [6,] 0.7031708 0.704707
 [7,] 0.8130203 0.814054
 [8,] 0.9118814 0.912410
 [9,] 1.0000000 1.000000
> 
> cbind(eecf(y)(z), round(ecMBBEFD(z, 3, 2), 6))
           [,1]     [,2]
 [1,] 0.0000000 0.000000
 [2,] 0.2086718 0.208350
 [3,] 0.3723734 0.371587
 [4,] 0.5086622 0.507853
 [5,] 0.6272755 0.626214
 [6,] 0.7327380 0.731859
 [7,] 0.8284908 0.828035
 [8,] 0.9170379 0.916909
 [9,] 1.0000000 1.000000
> 
> #test mean
> c(mean(x), mMBBEFD(1, 2, 1/5))
[1] 0.7584131 0.7590979
> 
> c(mean(y), mMBBEFD(1, 3, 2))
[1] 0.5166914 0.5169925
> 
> c(mean(x3), mMBBEFD(1, 3, 1))
[1] 0.5499008 0.5493061
> c(mean(x4), mMBBEFD(1, 3, 1/3))
[1] 0.6081893 0.6068262
> c(mean(x5), mMBBEFD(1, 3, 0))
[1] 1 1
> 
> #test E(X^2)
> c(mean(x^2), mMBBEFD(2, 2, 1/5))
[1] 0.6746952 0.6741345
> 
> c(mean(y^2), mMBBEFD(2, 3, 2))
[1] 0.4251654 0.4259820
> 
> c(mean(x3^2), mMBBEFD(2, 3, 1))
[1] 0.4503140 0.4506939
> c(mean(x4^2), mMBBEFD(2, 3, 1/3))
[1] 0.4994786 0.4978878
> c(mean(x5^2), mMBBEFD(2, 3, 0))
[1] 1 1
> 
> #total loss
> etl(x)
[1] 0.502
> tlMBBEFD(2, 1/2)
[1] 0.5
> 
> 
> etl(y)
[1] 0.332
> tlMBBEFD(3, 2)
[1] 0.3333333
> 
> 
> #test quantile
> cbind(quantile(x, probs=0:10/10), round(qMBBEFD(0:10/10, 2, 1/5), 6))
             [,1]     [,2]
0%   0.0007875513 0.000000
10%  0.2213643181 0.228480
20%  0.4217237630 0.430677
30%  0.6220458908 0.620421
40%  0.8123478423 0.807290
50%  1.0000000000 1.000000
60%  1.0000000000 1.000000
70%  1.0000000000 1.000000
80%  1.0000000000 1.000000
90%  1.0000000000 1.000000
100% 1.0000000000 1.000000
> cbind(quantile(y, probs=0:10/10), round(qMBBEFD(0:10/10, 3, 2), 6))
             [,1]     [,2]
0%   0.0001520552 0.000000
10%  0.0403069138 0.040642
20%  0.0962171587 0.093109
30%  0.1653841935 0.163499
40%  0.2613723876 0.263034
50%  0.4169472348 0.415037
60%  0.6730945360 0.678072
70%  1.0000000000 1.000000
80%  1.0000000000 1.000000
90%  1.0000000000 1.000000
100% 1.0000000000 1.000000
> 
> 
> z <- seq(0, 1, length=101)
> plot(z, pMBBEFD(z, 3, 2), type="l", ylim=c(0, 1-tlMBBEFD(3, 2)))
> 
> plot(z, qMBBEFD(z, 3, 2), type="l", xlim=c(0, 1-tlMBBEFD(3, 2)))
> 
> 
> 
> #test density
> 
> integrate(dMBBEFD, 0, 1, g=2, b=1/5)
0.5 with absolute error < 5.6e-15
> pMBBEFD(1-1e-6, 2, 1/5)
[1] 0.4999995
> 
> 
> integrate(dMBBEFD, 0, 1, g=3, b=2)
0.6666667 with absolute error < 2.9e-09
> pMBBEFD(1-1e-6, 3, 2)
[1] 0.6666665
> 
> 
> z <- 0:8/8
> 
> d <- function(z) approxfun(density(x)$x, density(x)$y)(z)
> cbind(d(z), dMBBEFD(z, 2, 1/5))
           [,1]      [,2]
 [1,] 0.2198903 0.4023595
 [2,] 0.4482357 0.4414640
 [3,] 0.4897915 0.4763761
 [4,] 0.4963786 0.5049045
 [5,] 0.4959501 0.5250597
 [6,] 0.5107150 0.5353281
 [7,] 0.5326246 0.5349045
 [8,] 0.6151932 0.5238226
 [9,] 4.7080922 0.5000000
> 
> d <- function(z) approxfun(density(y)$x, density(y)$y)(z)
> cbind(d(z), dMBBEFD(z, 3, 2))
           [,1]      [,2]
 [1,] 1.0835571 2.7725887
 [2,] 1.4615322 1.4330416
 [3,] 0.9107230 0.8706456
 [4,] 0.5882580 0.5826387
 [5,] 0.4346632 0.4157398
 [6,] 0.3220319 0.3104807
 [7,] 0.2378861 0.2398749
 [8,] 0.3907126 0.1902387
 [9,] 2.4204885 0.3333333
> 
> z <- sort(c(1, seq(-0.1,1.1, length=101)))
> plot(density(x), ylim=c(0,1))
> lines(z, dMBBEFD(z, 2, 1/5), col="red")
> 
> 
> plot(density(y), ylim=c(0,3))
> lines(z, dMBBEFD(z, 3, 2), col="red")
> 
> proc.time()
   user  system elapsed 
   1.04    0.17    1.20