simpleqt opened a new pull request, #140:
URL: https://github.com/apache/commons-statistics/pull/140

   `CauchyDistribution.inverseCumulativeProbability(p)` computes `location + 
scale * Math.tan(Math.PI * (p - 0.5))` (`CauchyDistribution.java:139`, and the 
survival function at `:156`). For a probability close to 0 or 1 the subtraction 
of one half cancels the significant digits of `p`: the product `Math.PI * (p - 
0.5)` is rounded to ± π/2, so the result saturates at ± 1.633 
&times; 10<sup>16</sup> &gamma; whatever the probability is:
   
   ```java
   final CauchyDistribution dist = CauchyDistribution.of(0.5, 2);
   dist.inverseCumulativeProbability(1e-20);           // 
-3.266247870639074e16, expected -6.366197723675814e19
   dist.inverseCumulativeProbability(1 - 1e-16);       //  
3.957875932190438e15, expected  5.73416113922266e15
   dist.inverseCumulativeProbability(1e-300);          // 
-3.266247870639074e16, expected -6.366197723675814e299
   ```
   
   This is not only a saturation of the extreme tail. The accuracy decreases 
from the middle of the tail: at `p = 1e-5` the value is out by 5e-12 relative 
(at `p = 1e-10` by 5e-7); the function is not monotonic, since consecutive 
probabilities below `1e-18` can return the identical value; and the median is 
offset by a rounding of &pi; (1.85e-16 relative for the scale 2.1, location 1.2 
distribution). A caller using the quantile function for a likelihood, a p-value 
or a rejection threshold in the tail gets a value that stops responding to the 
input and can decrease.
   
   Fix: use the equivalent forms valid for the two tails, `tan(pi (p - 1/2)) = 
-cot(pi p)` for `p < 1/2` and `cot(pi (1 - p))` for `p > 1/2`, which need no 
cancellation. The cotangent is evaluated as `1 / tan(x)` or `cos(x) / sin(x)` 
(the latter avoids the relative error of `sin(x)` for small `x`), and is 
exactly zero at the median. `inverseSurvivalProbability` shares the helper and 
is fixed by the same change.
   
   Validation:
   
   - `CauchyDistributionTest#testInverseProbabilityTails` checks 13 
probabilities from `1e-300` to `1 - 2^-53` against `scipy.stats.cauchy`, 
including both tails; on the pre-fix base it fails at `p = 1e-5` (`expected: 
<-63661.47721581417> but was: <-63661.47721549608>`).
   - `CauchyDistributionTest#testInverseProbabilityTailsAreMonotonic` checks 
strict monotonicity from `p = 0.1` down to `1e-300` and up to `1 - 1e-16`; on 
the pre-fix base it fails below `p = 1e-18`.
   - All 58 tests of `CauchyDistributionTest` pass, and so does the full 
reactor (`mvn test`: distribution 16247, descriptive 31001, inference 7434 
tests, 0 failures); checkstyle reports 0 violations.
   
   No `src/changes/changes.xml` entry is included because the action format 
requires a JIRA issue key and I do not have a JIRA account. Happy to add one, 
e.g.:
   
   ```xml
   <action dev="simpleqt" due-to="simpleqt" type="fix" 
issue="STATISTICS-xxx">CauchyDistribution: Fix the accuracy of the inverse 
probability functions in the tails, which saturated at 1.633e16 times the scale 
for probabilities close to 0 or 1.</action>
   ```
   
   AI assistance used to diagnose and fix this defect; the analysis and the 
change have been reviewed and verified as described above.
   


-- 
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.

To unsubscribe, e-mail: [email protected]

For queries about this service, please contact Infrastructure at:
[email protected]

Reply via email to