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 GB1 distribution
> 
> 
> 
> 
> #integral of the density
> integrate(dgbeta, 0, 1, shape0=1, shape1=3, shape2=3/2)
1 with absolute error < 3.3e-09
> integrate(dgbeta, 0, 1, shape0=1/2, shape1=3, shape2=3/2)
1 with absolute error < 3.2e-10
> integrate(dgbeta, 0, 1, shape0=2, shape1=3, shape2=3/2)
1 with absolute error < 2.4e-08
> 
> z <- 0:10/10
> cbind(pi*z^(3*pi-1)*(1-z^pi)^(3/2-1)/beta(3, 3/2), dgbeta(z, pi, 3, 3/2))
              [,1]         [,2]
 [1,] 0.000000e+00 0.000000e+00
 [2,] 7.749692e-08 7.749692e-08
 [3,] 2.655611e-05 2.655611e-05
 [4,] 8.018270e-04 8.018270e-04
 [5,] 8.894060e-03 8.894060e-03
 [6,] 5.649280e-02 5.649280e-02
 [7,] 2.491641e-01 2.491641e-01
 [8,] 8.384937e-01 8.384937e-01
 [9,] 2.233348e+00 2.233348e+00
[10,] 4.504939e+00 4.504939e+00
[11,] 0.000000e+00 0.000000e+00
> cbind(dbeta(z^pi, 3, 3/2)*pi*z^(pi-1), dgbeta(z, pi, 3, 3/2))
              [,1]         [,2]
 [1,] 0.000000e+00 0.000000e+00
 [2,] 7.749692e-08 7.749692e-08
 [3,] 2.655611e-05 2.655611e-05
 [4,] 8.018270e-04 8.018270e-04
 [5,] 8.894060e-03 8.894060e-03
 [6,] 5.649280e-02 5.649280e-02
 [7,] 2.491641e-01 2.491641e-01
 [8,] 8.384937e-01 8.384937e-01
 [9,] 2.233348e+00 2.233348e+00
[10,] 4.504939e+00 4.504939e+00
[11,] 0.000000e+00 0.000000e+00
> cbind(pbeta(z^pi, 3, 3/2), pgbeta(z, pi, 3, 3/2))
              [,1]         [,2]
 [1,] 0.000000e+00 0.000000e+00
 [2,] 8.223421e-10 8.223421e-10
 [3,] 5.639895e-07 5.639895e-07
 [4,] 2.559710e-05 2.559710e-05
 [5,] 3.802695e-04 3.802695e-04
 [6,] 3.044319e-03 3.044319e-03
 [7,] 1.634905e-02 1.634905e-02
 [8,] 6.588056e-02 6.588056e-02
 [9,] 2.110280e-01 2.110280e-01
[10,] 5.445815e-01 5.445815e-01
[11,] 1.000000e+00 1.000000e+00
> cbind(qbeta(z, 3, 3/2)^(1/pi), qgbeta(z, pi, 3, 3/2))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.7337871 0.7337871
 [3,] 0.7949537 0.7949537
 [4,] 0.8344269 0.8344269
 [5,] 0.8646673 0.8646673
 [6,] 0.8898449 0.8898449
 [7,] 0.9119622 0.9119622
 [8,] 0.9322443 0.9322443
 [9,] 0.9516895 0.9516895
[10,] 0.9716348 0.9716348
[11,] 1.0000000 1.0000000
> 
> 
> #log argument
> 
> cbind(dbeta(z, 3, 3/2, log=TRUE), dgbeta(z, 1, 3, 3/2, log=TRUE))
            [,1]       [,2]
 [1,]       -Inf        NaN
 [2,] -2.7764788 -2.7764788
 [3,] -1.4490760 -1.4490760
 [4,] -0.7049115 -0.7049115
 [5,] -0.2066226 -0.2066226
 [6,]  0.1485037  0.1485037
 [7,]  0.4015750  0.4015750
 [8,]  0.5660353  0.5660353
 [9,]  0.6303656  0.6303656
[10,]  0.5193581  0.5193581
[11,]       -Inf       -Inf
> cbind(log(dgbeta(z, pi, 3, 3/2, log=FALSE)), dgbeta(z, pi, 3, 3/2, log=TRUE))
             [,1]        [,2]
 [1,]        -Inf        -Inf
 [2,] -16.3730277 -16.3730277
 [3,] -10.5362506 -10.5362506
 [4,]  -7.1286177  -7.1286177
 [5,]  -4.7223717  -4.7223717
 [6,]  -2.8736422  -2.8736422
 [7,]  -1.3896434  -1.3896434
 [8,]  -0.1761482  -0.1761482
 [9,]   0.8035017   0.8035017
[10,]   1.5051743   1.5051743
[11,]        -Inf        -Inf
> cbind(log(pgbeta(z, pi, 3, 3/2, log=FALSE)), pgbeta(z, pi, 3, 3/2, log=TRUE))
             [,1]        [,2]
 [1,]        -Inf        -Inf
 [2,] -20.9188646 -20.9188646
 [3,] -14.3882302 -14.3882302
 [4,] -10.5730313 -10.5730313
 [5,]  -7.8746304  -7.8746304
 [6,]  -5.7944781  -5.7944781
 [7,]  -4.1135857  -4.1135857
 [8,]  -2.7199119  -2.7199119
 [9,]  -1.5557645  -1.5557645
[10,]  -0.6077377  -0.6077377
[11,]   0.0000000   0.0000000
> cbind(log(pgbeta(z, pi, 3, 3/2, log=FALSE, lower=TRUE)), pgbeta(z, pi, 3, 3/2, log=TRUE, lower=TRUE))
             [,1]        [,2]
 [1,]        -Inf        -Inf
 [2,] -20.9188646 -20.9188646
 [3,] -14.3882302 -14.3882302
 [4,] -10.5730313 -10.5730313
 [5,]  -7.8746304  -7.8746304
 [6,]  -5.7944781  -5.7944781
 [7,]  -4.1135857  -4.1135857
 [8,]  -2.7199119  -2.7199119
 [9,]  -1.5557645  -1.5557645
[10,]  -0.6077377  -0.6077377
[11,]   0.0000000   0.0000000
> cbind(qgbeta(z, pi, 3, 3/2, log=FALSE), qgbeta(log(z), pi, 3, 3/2, log=TRUE))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.7337871 0.7337871
 [3,] 0.7949537 0.7949537
 [4,] 0.8344269 0.8344269
 [5,] 0.8646673 0.8646673
 [6,] 0.8898449 0.8898449
 [7,] 0.9119622 0.9119622
 [8,] 0.9322443 0.9322443
 [9,] 0.9516895 0.9516895
[10,] 0.9716348 0.9716348
[11,] 1.0000000 1.0000000
> cbind(qgbeta(z, pi, 3, 3/2, log=FALSE, lower=TRUE), qgbeta(log(z), pi, 3, 3/2, log=TRUE, lower=TRUE))
           [,1]      [,2]
 [1,] 0.0000000 0.0000000
 [2,] 0.7337871 0.7337871
 [3,] 0.7949537 0.7949537
 [4,] 0.8344269 0.8344269
 [5,] 0.8646673 0.8646673
 [6,] 0.8898449 0.8898449
 [7,] 0.9119622 0.9119622
 [8,] 0.9322443 0.9322443
 [9,] 0.9516895 0.9516895
[10,] 0.9716348 0.9716348
[11,] 1.0000000 1.0000000
> 
> #RNG
> n <- 1e4
> x <- rgbeta(n, shape0=2, shape1=3, shape2=3/2)
> y <- rgbeta(n, shape0=pi, shape1=3, shape2=3/2)
> 
> #test density
> z <- 0:100/100
> 
> plot(density(x)); lines(z, dgbeta(z, shape0=2, shape1=3, shape2=3/2), col="red")
> plot(density(y)); lines(z, dgbeta(z, shape0=pi, shape1=3, shape2=3/2), col="red")
> 
> #mode
> modeGB1 <- function(shape0, shape1, shape2)
+ {
+   if(shape1+shape2-1/shape0>1)
+     ((shape1-1/shape0)/(shape1+shape2-1/shape0-1))^(1/shape0)
+   else
+     NaN
+ }
> c(modeGB1(2, 3, 3/2), density(x)$x[which.max(density(x)$y)])
[1] 0.9128709 0.9028143
> c(modeGB1(pi, 3, 3/2), density(y)$x[which.max(density(y)$y)])
[1] 0.9470343 0.9492062
> 
> #test CDF
> z <- 0:10/10
> cbind(ecdf(x)(z), pgbeta(z, shape0=2, shape1=3, shape2=3/2))
        [,1]         [,2]
 [1,] 0.0000 0.000000e+00
 [2,] 0.0000 2.179280e-06
 [3,] 0.0000 1.378829e-04
 [4,] 0.0013 1.539860e-03
 [5,] 0.0071 8.403941e-03
 [6,] 0.0301 3.079579e-02
 [7,] 0.0856 8.710400e-02
 [8,] 0.2035 2.041266e-01
 [9,] 0.4080 4.107520e-01
[10,] 0.7141 7.146727e-01
[11,] 1.0000 1.000000e+00
> 
> cbind(ecdf(y)(z), pgbeta(z, shape0=pi, shape1=3, shape2=3/2))
        [,1]         [,2]
 [1,] 0.0000 0.000000e+00
 [2,] 0.0000 8.223421e-10
 [3,] 0.0000 5.639895e-07
 [4,] 0.0000 2.559710e-05
 [5,] 0.0005 3.802695e-04
 [6,] 0.0030 3.044319e-03
 [7,] 0.0142 1.634905e-02
 [8,] 0.0647 6.588056e-02
 [9,] 0.2107 2.110280e-01
[10,] 0.5363 5.445815e-01
[11,] 1.0000 1.000000e+00
> 
> #plot(ecdf(x)); lines(z, pgbeta(z, shape0=2, shape1=3, shape2=3/2), col="red")
> #plot(ecdf(y)); lines(z, pgbeta(z, shape0=pi, shape1=3, shape2=3/2), col="red")
> 
> #mean
> c(mean(x), mgbeta(1, shape0=2, shape1=3, shape2=3/2))
[1] 0.8060811 0.8053399
> c(mean(y), mgbeta(1, shape0=pi, shape1=3, shape2=3/2))
[1] 0.8691161 0.8681286
> 
> #raw moment
> for(i in 2:4)
+ {
+   cat("E(X^", i, ")\n", sep="")
+   print(c(mean(x^i), mgbeta(i, shape0=2, shape1=3, shape2=3/2)))
+   print(c(mean(y^i), mgbeta(i, shape0=pi, shape1=3, shape2=3/2)))
+ }
E(X^2)
[1] 0.6675872 0.6666667
[1] 0.7645970 0.7629595
E(X^3)
[1] 0.5646389 0.5637379
[1] 0.6795740 0.6774769
E(X^4)
[1] 0.4856545 0.4848485
[1] 0.6093146 0.6068729
> 
> #test limited expected value
> d <- 1/2
> s0 <- 2
> s1 <- 3
> s2 <- 3/2
> 
> mean(pmin(x, d))
[1] 0.4978821
> 
> f <- function(x, d, shape0) dgbeta(x, shape0=shape0, shape1=3, shape2=3/2)*pmin(x,d)
> integrate(f, 0, 1, d=d, shape0=s0)
0.4977446 with absolute error < 5e-06
> 
> 
> 
> #test EC
> f <- function(x, d, shape0) dgbeta(x, shape0=shape0, shape1=3, shape2=3/2)*pmin(x, d)/mgbeta(1, shape0=shape0, shape1=3, shape2=3/2)
> F <- function(d, shape0) sapply(d, function(d) integrate(f, 0, 1, d=d, shape0=shape0)$value)
> 
> cbind(eecf(x)(z), ecgbeta(z, shape0=2, shape1=3, shape2=3/2), F(z, shape0=2))
           [,1]      [,2]      [,3]
 [1,] 0.0000000 0.0000000 0.0000000
 [2,] 0.1240570 0.1241711 0.1241712
 [3,] 0.2481140 0.2483374 0.2483373
 [4,] 0.3721111 0.3724309 0.3724291
 [5,] 0.4957368 0.4960794 0.4960793
 [6,] 0.6176576 0.6180553 0.6180553
 [7,] 0.7350465 0.7353826 0.7353816
 [8,] 0.8420269 0.8422583 0.8422585
 [9,] 0.9292422 0.9293007 0.9293013
[10,] 0.9845859 0.9844530 0.9844524
[11,] 1.0000000 1.0000000 1.0000000
> 
> cbind(eecf(y)(z), ecgbeta(z, shape0=pi, shape1=3, shape2=3/2), F(z, shape0=pi))
           [,1]      [,2]      [,3]
 [1,] 0.0000000 0.0000000 0.0000000
 [2,] 0.1150594 0.1151903 0.1151903
 [3,] 0.2301189 0.2303806 0.2303806
 [4,] 0.3451783 0.3455701 0.3455700
 [5,] 0.4602226 0.4607443 0.4607441
 [6,] 0.5751299 0.5757815 0.5757815
 [7,] 0.6893787 0.6900359 0.6900363
 [8,] 0.8004348 0.8010489 0.8010442
 [9,] 0.9009795 0.9016246 0.9016241
[10,] 0.9749285 0.9755137 0.9755141
[11,] 1.0000000 1.0000000 1.0000000
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.06    0.20    1.25