三次样条函数和插值函数是否都很小？

Question

我在scipy中使用三次样条函数来平滑一些概率密度函数（pdf），我遇到了一个问题：当y值很小时，scipy.interpolate.InterpolatedUnivariateSpline或{}没有通过所有的结

import math
import pandas as pd
import numpy as np
from scipy.stats import norm
from scipy.interpolate import InterpolatedUnivariateSpline, interp1d


def main():
    breadth = 0.60
    stepsize = 0.01
    mean = 0.0
    stdev = 0.02

    steps = 2. * breadth / stepsize + 1
    retRange = np.linspace(-breadth, breadth, num=steps)
    probs = {}
    for ret in retRange:
        retup = ret + stepsize / 2.
        retdown = ret - stepsize / 2.
        probup = norm.cdf(retup, loc=mean, scale=stdev)
        probdown = norm.cdf(retdown, loc=mean, scale=stdev)
        prob = probup - probdown
        probs[ret] = prob
    probs = pd.Series(probs, name='PDF')
    probs = probs / probs.sum()

    spl = InterpolatedUnivariateSpline(probs.index, probs, k=3, ext='zeros')
    tmp = pd.Series(index=retRange, data=retRange, name='SmoothedPDF').map(spl)

    out = pd.concat([probs, tmp], axis=1)
    out['diff'] = out['PDF'] - out['SmoothedPDF']
    print out

if __name__ == '__main__':
    main()

问题是输出数据帧显示了差异。我可以用interp1d代替InterpolateUnivariateSpline，得到同样的问题（尽管结果不同）。值得注意的是，在节点处插值的y值的差值可能大于原始y值。在

^{pr2}$

是的，我正在研究一些涉及非常低概率的事情，因此规模确实（有时）那么小。在

如果是精度问题，我可以做一些重新缩放，但在这样做之前，我想确认这是问题所在。在

作为旁注，使用线性插值（k=1）并没有显示出这个问题。在

Answer 1

问题似乎是你尝试使用的样条插值器全局地使用插值，这对你的输入非常不合适。由于指数截止，很多值由于双重精度而精确为零；在这个范围内，样条插值器会给您提供噪声。在

这是一个问题的数值表示，通过在半对数标度上绘制PDF及其样条插值（请注意，如果PDF为零，则在对数图中会缺少点）：

import matplotlib.pyplot as plt

plt.figure()
plt.semilogy(out.index, out.PDF, 'r-', out.index, out.SmoothedPDF.abs(), 'bo ',
             markerfacecolor='none') # note the abs for log
plt.legend(['PDF', 'SmoothedPDF (univariate spline)'])
plt.show()

首先，我需要对插值值调用.abs()，以便绘制出负的值（我们插值了一个非零函数！）。其次，正如我已经指出的，x=0.2以上的范围应该是缺失的，因为函数值在这里正好是零。相反，我们有1e-50阶的噪声，就像在小x极限上一样。总之，我们能看到的是，如果你的数据小于1e-30，你就没有机会了。在

我的第一个想法是使用griddata，但是它显示了同样糟糕的结果（即使在这种情况下，数值误差要小得多）。那么使用一个局部插值器怎么样？输入piecewise cubic Hermite spline。为了演示插值和数值稳定性：

^{pr2}$

稳定性证明：

>>> out['diff'].any()
False

误差标度对数图：

如你所见，协议几乎是完美的。请注意，在原始数据正好为0的情况下，明显缺少值。在小x范围内，什么样的噪声保持在1e-100以下，你必须学会如何生活。在

插值结果的实际绘图：