CODEWITHSUNDEEP
pythonfunctionnumpyscipy
Mc Missile
Mc Missile

Reputation: 735

Invert a function in python without solving for it

I have a function that calculates the cumulative distributive function(cdf) of a skewed gaussian as follows.

from scipy import specials

def cdf(x,location,scale,a):
    temp = (x-location)/abs(scale)
    return 0.5*(1+special.erf(temp/np.sqrt(2))) - 2.0*special.owens_t(temp,a)

How do I find the value of x(accurately) that gives a cdf of 99.9987. My current method incorporates using linspace between two limits, calculating the cdf for the linspace array and finding the value that is closest to 99.9987. However, this method is not very robust and breaks easily due to cdf not being a linear function.

Upvotes: 1

Views: 1983

Answers (1)

xdze2
xdze2

Reputation: 4151

You could use one of the scalar root finding solver proposed by scipy (see for example the bisection method):

import numpy as np
from scipy.optimize import root_scalar
from scipy import special

def cdf(x, location, scale, a):
    temp = (x - location)/abs(scale)
    return 0.5*(1 + special.erf(temp/np.sqrt(2))) - 2.0*special.owens_t(temp, a)

# Root finding
location, scale, a = 50, 2, 5
target = 0.999987
sol = root_scalar(lambda x, *args: cdf(x, *args) - target,
                  bracket=(-6*scale + location, 6*scale + location),
                  args=(location, scale, a))

print(sol.flag)

# Graph
x = np.linspace(-6*scale + location, 6*scale + location, 123)
y = cdf(x, location, scale, a)

plt.plot(x, y);
plt.plot(sol.root, target, '.r');
plt.xlabel('x'); plt.ylabel('cdf');

which gives:

resulting graph

Upvotes: 1

Related Questions