search

swihart / mvgb

Multivariate probabilities via the Genz and Bretz method
GNU Lesser General Public License v2.1
0 stars 0 forks source link

loss of monotonicity and positivity in libstable #2

Open swihart opened 2 years ago

swihart commented 2 years ago

While for this package it has been affecting the code in F77.c, the loss of monotonicity and positivity can be demonstrated with following R code:


set.seed(10)
n <- 250000
randunif <- sort(runif(n))
the_draw <- libstableR::stable_q(randunif, 
                                 c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                                 parametrization = 1L)

randunif[the_draw < 0]
the_draw[the_draw < 0]

## smallest/biggest of the random uniform
head(sort(randunif))
tail(sort(randunif))

## smallest/biggest of the draw
head(the_draw)
tail(the_draw)

## sqrt() those
head(sqrt(the_draw))
tail(sqrt(the_draw))

## now reciprocate
head(1/sqrt(the_draw))
tail(1/sqrt(the_draw))

## goes negative if tol too tight?
set.seed(10)
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-6))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-7))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-8))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-9))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-10))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-11))
libstableR::stable_q(0.9997662,
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=c(1e-12))

## if tol is too tight we lose monotonicity?
libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-6)
libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-7)
libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-8)
libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-9)
libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-10)

stabledist::qstable(seq(0.9998,0.9999,0.00001), 
                    1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0,
                    pm = 1L,
                    tol=1e-12)

libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
                     c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
                     parametrization = 1L, tol=1e-12)

library(stable, lib.loc="~/RLIBS/")
stable::qstable(seq(0.9998,0.9999,0.00001), 
                1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0,
                param = 1
)

which gives:

> set.seed(10)
> n <- 250000
> randunif <- sort(runif(n))
> the_draw <- libstableR::stable_q(randunif, 
+                                  c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                                  parametrization = 1L)
> 
> randunif[the_draw < 0]
[1] 0.9997662
> the_draw[the_draw < 0]
[1] -116.5775
> 
> ## smallest/biggest of the random uniform
> head(sort(randunif))
[1] 3.079884e-06 5.020527e-06 1.001381e-05 1.006923e-05 1.136865e-05 1.339754e-05
> tail(sort(randunif))
[1] 0.9999678 0.9999724 0.9999778 0.9999792 0.9999795 0.9999961
> 
> ## smallest/biggest of the draw
> head(the_draw)
[1] 0.8017419 0.8079077 0.8174805 0.8175604 0.8193301 0.8217647
> tail(the_draw)
[1]  17578.95  21020.87  27148.07  29367.69  29817.23 207702.20
> 
> ## sqrt() those
> head(sqrt(the_draw))
[1] 0.8954004 0.8988368 0.9041463 0.9041905 0.9051686 0.9065124
Warning message:
In sqrt(the_draw) : NaNs produced
> tail(sqrt(the_draw))
[1] 132.5856 144.9858 164.7667 171.3700 172.6767 455.7436
Warning message:
In sqrt(the_draw) : NaNs produced
> 
> ## now reciprocate
> head(1/sqrt(the_draw))
[1] 1.116819 1.112549 1.106016 1.105962 1.104767 1.103129
Warning message:
In sqrt(the_draw) : NaNs produced
> tail(1/sqrt(the_draw))
[1] 0.007542296 0.006897228 0.006069187 0.005835326 0.005791170 0.002194216
Warning message:
In sqrt(the_draw) : NaNs produced
> 
> 
> 
> ## goes negative if tol too tight?
> set.seed(10)
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-6))
[1] 1707.877
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-7))
[1] 2046.2
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-8))
[1] 2046.2
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-9))
[1] -116.5717
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-10))
[1] -116.5717
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-11))
[1] -116.5716
> libstableR::stable_q(0.9997662,
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=c(1e-12))
[1] -116.5717
> 
> 
> ## if tol is too tight we lose monotonicity?
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-6)
 [1] 2051.263 2178.611 2321.454 2482.702 2666.022 2876.109 3119.081 3403.007 3738.807
[10] 4141.576 4632.742
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-7)
 [1] 2051.263 2178.611 2321.454 2482.702 2666.022 2876.109 3119.081 3403.007 3738.807
[10] 4141.576 4632.742
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-8)
 [1] 2051.263 2178.611 2321.454 2482.702 2666.022 2876.109 3119.081 3403.007 3738.807
[10] 4141.576 4632.742
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-9)
 [1] 3312.472 1143.143 4932.699 4046.233 3974.338 5010.710 3738.837 4078.938 4481.211
[10] 4143.174 4633.845
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-10)
 [1] 3312.472 1143.143 4932.699 4046.233 3974.338 5010.710 3738.837 4078.938 4481.211
[10] 4143.174 4633.845
> 
> 
> 
> 
> stabledist::qstable(seq(0.9998,0.9999,0.00001), 
+                     1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0,
+                     pm = 1L,
+                     tol=1e-12)
 [1] 4395.411 4395.411 4395.411 4395.411 4395.411 4395.411 4395.411 4395.411 4395.411
[10] 4395.411 4395.411
> 
> 
> libstableR::stable_q(seq(0.9998,0.9999,0.00001), 
+                      c(1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0), 
+                      parametrization = 1L, tol=1e-12)
 [1] 3312.472 1143.143 4932.699 4046.233 3974.338 5010.710 3738.837 4078.938 4481.211
[10] 4143.174 4633.845
> 
> 
> library(stable, lib.loc="~/RLIBS/")
> stable::qstable(seq(0.9998,0.9999,0.00001), 
+                 1.7/2, 1, 2*cos(1.7/2*pi/2)^(2/1.7),0,
+                 param = 1
+ )
 [1]  5239.963  5565.609  5930.815  6343.012  6811.576  7348.513  7969.429  8694.977
 [9]  9553.047 10582.200 11837.231
swihart commented 2 years ago

Seems like Robust Analysis values are higher than the other CRAN alternatives.

My workaround for now is to basically draw random variables in F77.c (as opposed to stable_q(uniform)) and crank up the MAXPTS to make the returned inform a value of 0. Preliminary testing shows that despite the increase in MAXPTS this approach may be 50% or so faster.