autograd.scipy.stats.norm.ppf not implemented

[Johneinsteinwong]

hi,

i would like to get the grad of the following function wrt to the 1st argument :

def Phi(s0,z,b):
  phi = norm.sf( norm.ppf(s0) + np.dot(z,b) )
  return phi

small_phi = egrad(Phi,0)

However, i found that norm.ppf has not been implemented in autograd.scipy.stats. But it is available in scipy:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html

Inside the autograd.scipy.stats.norm.py script, i found the wrappers for cdf and others function like pdf:

defvjp(cdf,
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(x, lambda g: g * pdf(x, loc, scale)) ,
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(loc, lambda g: -g * pdf(x, loc, scale)),
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(scale, lambda g: -g * pdf(x, loc, scale)*(x-loc)/scale))

Could someone pls explain to me how does it work? How can i implement the wrapper for norm.ppf and norm.isf? Thanks

[CamDavidsonPilon]

Using your code snippet as an template, the structure looks like:

defvjp(cdf,
         <how to handle derivative w.r.t. argument 0, x>, 
         <how to handle derivative w.r.t. argument 1, location>,
         <how to handle derivative w.r.t. argument 2, scale>,

Most of the derivative code is boilerplate. Focusing on <how to handle derivative w.r.t. argument 0, x>, we can see, in code, that CDF'(x) = pdf(x):

lambda g: g * pdf(x, loc, scale)

(The g is there because of the chain rule)

Let's think about how ppf might look like that. ppf(q) = CDF^{-1}(q). From calculus, we know that (f^{-1})'(x) = 1/f'(f^{-1})(x). For ppf then,

defvjp(ppf,
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(x, lambda g: g * 1/pdf(ppf(x, loc, scale), loc, scale)) ,

That should work. For the other arguments, it will get more complicated, but the same idea applies. The below code is everything you need, and is the start of a PR to autograd library.

from scipy.stats import norm
from autograd.scipy.stats import norm as ag_norm
from autograd.extend import primitive, defvjp
from autograd.numpy.numpy_vjps import unbroadcast_f

ppf = primitive(norm.ppf)

defvjp(ppf,
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(x, lambda g: g * 1/ag_norm.pdf(ppf(x, loc, scale), loc, scale)) 
)

from autograd import grad

print(grad(ppf)(0.5))
print(grad(ppf)(0.0))

Does that give you enough information to do the rest?

[MauroCE]

@CamDavidsonPilon when I run your code I get RuntimeWarning: divide by zero encountered in double_scalars unbroadcast_f(x, lambda g: g * 1/ag_norm.pdf(ppf(x, loc, scale), loc, scale)). Is this normal?

[MauroCE]

I am trying to make this work for a function like this:

from numpy import exp
from scipy.stats import norm
from autograd.scipy.stats import norm as ag_norm
from autograd.extend import primitive, defvjp
from autograd.numpy.numpy_vjps import unbroadcast_f

ppf = primitive(norm.ppf)

defvjp(ppf,
       lambda ans, x, loc=0.0, scale=1.0:
       unbroadcast_f(x, lambda g: g * 1/ag_norm.pdf(ppf(x, loc, scale), loc, scale)) 
)

def f(thetau):
    a, b, g, k, *u = thetau
    z = ppf(u)
    return a + b*(1 + 0.8 * (1 - exp(-g * z)) / (1+exp(-g * z))) * ((1 - z)**k) * z

Here thetau is the concatenation between a vector of fixed length theta = [a, b, g, k] and a vector with variable length generated as u = rand(N) for some N. Thankfully the scipy ppf function broadcasts to a vector, but this doesn't seem to work here. My guess is that this is due to the use of unbroadcast_f. I tried removing it but nothing changes, I get this error:

TypeError: ufunc 'ndtri' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

[MauroCE]

Aaah I've solved it! The problem is that it doesn't work for a list, which is what *u gives. Therefore just adding u = np.array(u) fixes the problem.