1
$\begingroup$

I would like to know if there exists a package or how one can fast calculate the quantiles of a function within python, where the inverse of the function for calculating the quantile depends on the first and second derivative. More concretely.

Suppose we are given a $C^2([a, b])$ function $f$ and define the curvature $\kappa$ as

$$ \kappa(x) = \frac{|f''(x)|}{|1+(f'(x))^2|^{\frac{3}{2}}}$$

This is a positive function. I then normalize $\kappa$ s.t. $\int_a^b \kappa = 1$. For simplicity let's denote the normalized function with $g(x):=\frac{\kappa(x)}{\int_a^b \kappa}$.

Now I define the cumulative distribution function (CDF) $F$ as:

$$ F(x) = \int_a^x g(t) dt $$

we know that $F(a) = 0$ and $F(b) = 1$. The problem I want to solve is:

For given value $y$ of $F$ find $x$ s.t. $F(x) = y$, i.e. the inverse problem of $F$.

I'm interested in a solution within python (or any other similar language like Matlab, R). In my problem I only have data points of $f$ and need to approximate it. The following code snippet is an example for the $sin^2(x)$ curve.

import timeit
import numpy as np
from scipy.interpolate import UnivariateSpline
import scipy.integrate as integrate
from scipy.optimize import brentq


x = np.linspace(0, 6, 100)
interp = UnivariateSpline(x, np.sin(x)**2, k=4, s=0)

# calculate the normalization constant


def kappa(x):
    return abs(interp.derivative(2)(x))/abs((1+interp.derivative(1)(x)**2)**(1.5))


k = integrate.quad(kappa, 0, 6)[0]


def g(x):
    return kappa(x)/k


x = np.linspace(0, 6, 100)
interp2 = UnivariateSpline(x, g(x), k=4, s=0)


def F(x):
    return integrate.quad(interp2, 0, x)[0]

# get 20 test points
# not 0 and 1 are clear so we don't need
# to estimate these points.


N = 20
x2 = [(i - 1)/(N-1) for i in range(2, N)]

# first approach would be
start_time = timeit.default_timer()
y = [brentq(lambda u: F(u) - i, 0, 6) for i in x2]
elapsed = timeit.default_timer() - start_time
print(elapsed)

# just add the two default points

# y = [0] + y
# y = y.insert(len(y), 6)


def F2(x):
    return integrate.quad(g, 0, x)[0]


start_time = timeit.default_timer()
y = [brentq(lambda u: F2(u) - i, 0, 6) for i in x2]
elapsed = timeit.default_timer() - start_time
print(elapsed)

Note, I've calculated the inverse with two different methods. First I use another interpolation to approximate $\kappa$. In the second version I run it via directly using $\kappa$ where it has to evaluate the first and second derivative. If I'm running this on my pc I get the elapsed times:

0.449614048004
7.88159608841

What a huge different! Why is the second one so much faster? Now given what I want to do what is the best way to speed this up (even more than in case 1).

EDIT I've actually had a bug in the code. Since I've used python 2.7 I needed to import the division module from future. With this I get the following. Time is still a huge difference and values are very close for both system:

Time and values for first method:

2.38353395462
[0.1437681151139289, 0.35456065909107926, 1.1539903423627926, 1.3971602056486263, 1.5460856607363669, 1.686478225698852, 1.8751150788777717, 2.6398043923321652, 2.9355086219617466, 3.091753998559681, 3.230312179133905, 3.4029515812391478, 4.058035698689046, 4.470426379024948, 4.636764697040471, 4.774960068200856, 4.93583175345252, 5.273899782448736]

Time and values for second method

57.0844209194
[0.1437734552194811, 0.354572374885639, 1.1530574841452361, 1.3967524608790778, 1.545758378815468, 1.6861161414005308, 1.8745107127405394, 2.637806584725155, 2.9348568464428215, 3.0912665934557833, 3.229804581185758, 3.4021836211797964, 4.0530805425833485, 4.469683619221007, 4.636255181639525, 4.774458238739285, 4.935134505525859, 5.271168509499432]
$\endgroup$
5
  • 1
    $\begingroup$ @TommiBrander thx for your comment I switched from inverse problem to inverse. it's not a pure python question since it deals with finding the inverse. Hope that fits it better $\endgroup$
    – math
    Nov 10 '17 at 18:07
  • $\begingroup$ are the resulting answers for the two about the same? Can you show their associated inverse values? $\endgroup$
    – spektr
    Nov 10 '17 at 19:34
  • $\begingroup$ @C.Howard added the output. There was actually a small bug in my code (it used always $0$ for $x$). The finding is the same. Method 2 is much much slower. Values are very similar. Do you know how I can tune this? Even method 1 look extremely slow. $\endgroup$
    – math
    Nov 10 '17 at 19:51
  • $\begingroup$ Your post does not have to reflect its history. If anybody really needs to know, it’s publicly available anyway. Please edit your question to only contain what is relevant now, i.e., how it would be if you wrote it from scratch. In particular please fix all bugs in your example code, so everybody can easily reproduce your results. $\endgroup$
    – Wrzlprmft
    Nov 10 '17 at 20:52
  • 1
    $\begingroup$ I think your second one, $F_2(\cdot)$, is slower due to the fact you are computing the derivatives of the $\sin^2(x)$ spline over and over. The derivative computation seems to be the main major computation and so if you're seeing major run-time differences, the problem probably resides there. You should try to find a way to profile your code and see if you can identify what the most time consuming part is for $F_2(\cdot)$. $\endgroup$
    – spektr
    Nov 11 '17 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.