Percent point function (inverse of `CDF`) in pymc

Questions v5
,
1

I would like to calculate the inverse CDF of a Variable in pymc. I would typically use norm.ppf from Scipy, however within this function there are boolean operations and so I get the type error “Variables do not support boolean operations.”

I would normally get around this error by using aesara.ifelse, however here I am dealing with boolean operations within the package norm.ppf (which I believe is written in C or fortran).

Any ideas?

p.s., error can be replicated by simply running:

from scipy.stats import norm
import aesara

print(norm.ppf(0.5)) ## desired use of ppf (but with variable instead of float)
print(norm.ppf(aesara.tensor.as_tensor_variable(0.5))) ## fails

with error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 norm.ppf(aesara.tensor.as_tensor_variable(0.5))

File c:\Users\spearing\AppData\Local\Continuum\anaconda3\envs\pymc4_env\lib\site-packages\scipy\stats\_distn_infrastructure.py:2156, in rv_continuous.ppf(self, q, *args, **kwds)
   2154 _a, _b = self._get_support(*args)
   2155 cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
-> 2156 cond1 = (0 < q) & (q < 1)
   2157 cond2 = cond0 & (q == 0)
   2158 cond3 = cond0 & (q == 1)

File c:\Users\spearing\AppData\Local\Continuum\anaconda3\envs\pymc4_env\lib\site-packages\aesara\tensor\var.py:71, in _tensor_py_operators.__bool__(self)
     69     return True
     70 else:
---> 71     raise TypeError("Variables do not support boolean operations.")

TypeError: Variables do not support boolean operations.
1 Like
2

I think there is one, but I forgot. But you can use the pm.math.erf to compute it:

    def cdf(x):
        return 0.5 + 0.5 * pm.math.erf(x / pm.math.sqrt(2))

    def invcdf(x):
        return -pm.math.sqrt(2) * pm.math.erfinv(2*x)
1 Like
3

Thanks a lot!
The the cdf function works, but I think the invcdf function needs a little shift, but this works for me now.

def cdf(x):
    return 0.5 + 0.5 * pm.math.erf(x / pm.math.sqrt(2))

def invcdf(x):
    x = 0.5-x
    return -pm.math.sqrt(2) * pm.math.erfinv(2*x)
2 Likes
4

I’m glad it was helpful. By the way, I’ve finally found the functions in PyMC. they are:

inverse cdf: pm.math.probit 
cdf: pm.math.invprobit

You can see the source here: pymc.math — PyMC 4.1.7 documentation

Note that pm.math.probit uses the formula I posted before, without the shift you added. I don’t quite get why have you added that shift, the only different version of the formula I’ve seen (as in Wikipedia) is:

def invcdf(x):
        return pm.math.sqrt(2) * pm.math.erfinv(2*x - 1)

I have only tried the version I posted above, the one used by PyMC,and it usually works fine for me (I don’t know your specific case though).

1 Like
5

Ahh that’s even better, thanks!

I just added the shift so it lines up with the (0,1) scale for a standard Uniform variable:
image

That other version you mention is the same calculation but the transformation is inside the erfinv function, whereas I just transformed first.

6

I see, thanks for the clarification :slightly_smiling_face: