R Under development (unstable) (2024-10-16 r87241 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.12
Date:     
BugReport: https://github.com/spedygiorgio/mbbefd/issues


> 
> #test of GB1 distribution
> 
> #integral of the improper density
> integrate(doigbeta, 0, 1, shape0=1, shape1=3, shape2=3/2, p1=1/3)
0.6666667 with absolute error < 2.2e-09
> integrate(doigbeta, 0, 1, shape0=1/2, shape1=3, shape2=3/2, p1=1/3)
0.6666668 with absolute error < 0.00011
> integrate(doigbeta, 0, 1, shape0=2, shape1=3, shape2=3/2, p1=2/3)
0.3333333 with absolute error < 7.9e-09
> 
> 
> #RNG
> n <- 1e4
> x <- roigbeta(n, shape0=2, shape1=3, shape2=3/2, p1=1/3)
> y <- roigbeta(n, shape0=pi, shape1=3, shape2=3/2, p1=2/3)
> 
> c(etl(x), tloigbeta(shape0=2, shape1=3, shape2=3/2, p1=1/3))
[1] 0.3333000 0.3333333
> c(etl(y), tloigbeta(shape0=pi, shape1=3, shape2=3/2, p1=2/3))
[1] 0.6644000 0.6666667
> 
> #test CDF
> z <- 0:10/10
> cbind(ecdf(x)(z), poigbeta(z, shape0=2, shape1=3, shape2=3/2, p1=1/3))
        [,1]         [,2]
 [1,] 0.0000 0.000000e+00
 [2,] 0.0000 1.452854e-06
 [3,] 0.0002 9.192194e-05
 [4,] 0.0010 1.026573e-03
 [5,] 0.0056 5.602627e-03
 [6,] 0.0202 2.053053e-02
 [7,] 0.0596 5.806933e-02
 [8,] 0.1383 1.360844e-01
 [9,] 0.2764 2.738347e-01
[10,] 0.4752 4.764485e-01
[11,] 1.0000 1.000000e+00
> cbind(ecdf(y)(z), poigbeta(z, shape0=pi, shape1=3, shape2=3/2, p1=2/3))
        [,1]         [,2]
 [1,] 0.0000 0.000000e+00
 [2,] 0.0000 2.741140e-10
 [3,] 0.0000 1.879965e-07
 [4,] 0.0000 8.532368e-06
 [5,] 0.0000 1.267565e-04
 [6,] 0.0014 1.014773e-03
 [7,] 0.0058 5.449682e-03
 [8,] 0.0225 2.196019e-02
 [9,] 0.0717 7.034267e-02
[10,] 0.1848 1.815272e-01
[11,] 1.0000 1.000000e+00
> 
> 
> #mean
> c(mean(x), moigbeta(1, shape0=2, shape1=3, shape2=3/2, p1=1/3))
[1] 0.8698843 0.8702266
> c(mean(y), moigbeta(1, shape0=pi, shape1=3, shape2=3/2, p1=2/3))
[1] 0.9553983 0.9560429
> 
> #raw moment
> for(i in 2:4)
+ {
+   cat("E(X^", i, ")\n", sep="")
+   print(c(mean(x^i), moigbeta(i, shape0=2, shape1=3, shape2=3/2, p1=1/3)))
+   print(c(mean(y^i), moigbeta(i, shape0=pi, shape1=3, shape2=3/2, p1=2/3)))
+ }
E(X^2)
[1] 0.7773201 0.7777778
[1] 0.9198620 0.9209865
E(X^3)
[1] 0.7087100 0.7091586
[1] 0.8910033 0.8924923
E(X^4)
[1] 0.6561889 0.6565657
[1] 0.8671878 0.8689576
> 
> 
> #test EC
> cbind(eecf(x)(z), ecoigbeta(z, shape0=2, shape1=3, shape2=3/2, p1=1/3))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.1149578 0.1149126
 [3,] 0.2299134 0.2298222
 [4,] 0.3448175 0.3446868
 [5,] 0.4594474 0.4592770
 [6,] 0.5730509 0.5728352
 [7,] 0.6838473 0.6835254
 [8,] 0.7878653 0.7877674
 [9,] 0.8796336 0.8797731
[10,] 0.9520428 0.9521039
[11,] 1.0000000 1.0000000
> 
> cbind(eecf(y)(z), ecoigbeta(z, shape0=pi, shape1=3, shape2=3/2, p1=2/3))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.1046684 0.1045978
 [3,] 0.2093368 0.2091956
 [4,] 0.3140052 0.3137932
 [5,] 0.4186735 0.4183862
 [6,] 0.5232852 0.5229376
 [7,] 0.6276446 0.6272522
 [8,] 0.7310675 0.7305856
 [9,] 0.8312000 0.8307599
[10,] 0.9231080 0.9228566
[11,] 1.0000000 1.0000000
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.35    0.29    1.59