计算科学 - 如何快速估计用于计算分位数的导数 - 吾爱随笔录

我想知道是否存在一个包或者如何快速计算python中函数的分位数，其中用于计算分位数的函数的逆取决于一阶和二阶导数。更具体地说。

假设我们有一个 $C^2([a, b])$ 功能 $f$ 并定义曲率 $\kappa$ 作为

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

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

这是一个积极的功能。然后我规范化 st。为简单起见，让我们用来表示归一化函数。 $\kappa$ $\int_a^b \kappa = 1$ $g(x):=\frac{\kappa(x)}{\int_a^b \kappa}$

现在我将累积分布函数 (CDF)定义为： $F$

F (x) = \int_{a}^{x} g (t) d t

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

我们知道和。我要解决的问题是： $F(a) = 0$ $F(b) = 1$

对于的给定值找到 st，即的逆问题。 $y$ $F$ $x$ $F(x) = y$ $F$

我对 python（或任何其他类似的语言，如 Matlab、R）中的解决方案感兴趣。在我的问题中，我只有的数据点并且需要对其进行近似。以下代码片段是曲线的示例。 $f$ $sin^2(x)$

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)

注意，我用两种不同的方法计算了逆。首先，我使用另一个插值来近似。在第二个版本中，我直接使用运行它，它必须评估一阶和二阶导数。如果我在我的电脑上运行它，我会得到经过的时间： $\kappa$ $\kappa$

0.449614048004
7.88159608841

多么大的不同！为什么第二个这么快？现在考虑到我想做的事情，加快速度的最佳方法是什么（甚至比案例 1 还要多）。

编辑我实际上在代码中有一个错误。由于我使用了 python 2.7，我需要从未来导入除法模块。有了这个，我得到以下信息。时间仍然是一个巨大的差异，两个系统的值非常接近：

第一种方法的时间和价值：

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]

第二种方法的时间和价值

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]