BUG: ICDF implementation for the discrete geometric distribution fails some tests.

[gokuld]

Describe the issue:

The ICDF function for the discrete geometric distribution involves the ceiling function applied on a ratio of two values which is a floating point number. If this ratio is close to an integer, a small noise added to it even in the smallest digit of precision can sway the result of the ceiling function to the other side resulting in a wrong ICDF value.

The example code below is for an example of a Geometric distribution with p=0.99, and CDF value 0.9999 whose ICDF is expected to be 2. Within the ICDF function, defining the ratio a = pt.log1p(-value) / pt.log1p(-p)

    Ratio a before ceil: 2.000000000000025
    Ratio a after ceil: 3.0
    Expected value after ceil: 2

The tiny numerical error 000000000000025 contributes to the wrong value after ceil.

I also have a solution to this, which involves truncating the floating point value to a certain number of digits before using the ceiling function on it. This gets rid of the sensitivity to the smallest of noise. I will push a commit with this fix shortly.

Reproduceable code example:

import pytensor as pt
import pymc as pm
from pymc.logprob.basic import icdf
from pymc.pytensorf import inputvars
import numpy as np

value = 0.9999 # CDF value for n=2
p = 0.99 # Parameter p of the discrete Geometric distribution.

# The below ratio 'a' is used inside the ICDF function:
a = np.log1p(-value) / np.log1p(-p)

print(
    f"""
    Ratio before ceil: {a}
    Ratio after ceil: {np.ceil(a)}
    Expected value after ceil: 2
    """
)

dist = pm.Geometric.dist(p=p)

dist_icdf = icdf(dist, pt.tensor.type.TensorType(dtype='float64', shape=[])('value'))
dist_icdf_fn = pt.function(list(inputvars(dist_icdf)), dist_icdf)

assert dist_icdf_fn(value) == 2

Error message:

AssertionError

PyMC version information:

pytensor: 2.10.1
Python: 3.10.10
PyMC: latest commit 2a324bc

Context for the issue:

No response

[ricardoV94]

Floating point is tricky. Nice investigation.

Instead of truncating the number, maybe we can do as Scipy does? They manually check if if the logcdf disagrees, and tweak the value down:

https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/stats/_discrete_distns.py#L437-L440

[gokuld]

Created a pull request with a fix here #6671.

[gokuld]

Floating point is tricky. Nice investigation.

Yes, it is. Thank you. 👍

Instead of truncating the number, maybe we can do as Scipy does? They manually check if if the logcdf disagrees, and tweak the value down:

https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/stats/_discrete_distns.py#L437-L440

Let me get back after I check this.

[gokuld]

Yes the scipy approach sounds like a good idea, more computation (CDF is computed again) but also less error prone. Let me implement this and get back.

Is there any deadline for any release due to which this bug will be a blocker?

[ricardoV94]

No. I wouldn't count this as a blocker bug, or even a bug bug. More of a precision issue.

We will merge the solution whenever it's ready and the subsequent release will incorporate it. No rush

[gokuld]

I have pushed code with the suggested changes in #6671.