Кажется, проблема в том, что сплайн-интерполятор, который вы пытаетесь использовать, интерполирует глобально, что очень плохо подходит для вашего ввода.Из-за экспоненциального отсечения многие ваши значения равны нулю благодаря двойной точности;в этом диапазоне сплайн-интерполяторы дают вам шум.
Вот числовое представление проблемы путем построения вашего 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, точно так же как на пределе small-x.В целом, мы можем видеть, что если ваши данные меньше 1e-30, у вас нет шансов.
Моя первая идея состояла в том, чтобы использовать griddata, но он показывает такие же плохие результаты (хотячисловая ошибка намного меньше в этом случае).Так как насчет использования локального интерполятора?Введите кусочно кубический сплайн Эрмита .С немного измененной версией вашего кода, чтобы продемонстрировать как интерполяцию, так и числовую стабильность:
import math
import pandas as pd
import numpy as np
from scipy.stats import norm
from scipy.interpolate import pchip_interpolate
import matplotlib.pyplot as plt
def main():
breadth = 0.60
stepsize = 0.01
mean = 0.0
stdev = 0.02
steps = 2. * breadth / stepsize + 1
bigRange = np.linspace(-breadth, breadth, num=5*steps)
retRange = bigRange[::5]
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()
tmp_big = pd.Series(index=bigRange,
data=pchip_interpolate(probs.index.values, probs.values, bigRange),
name='SmoothedPDF_big') # denser interpolated case
tmp_check = pd.Series(index=retRange,
data=pchip_interpolate(probs.index.values, probs.values, retRange),
name='SmoothedPDF') # base points for checking stability
out = pd.concat([probs, tmp_check], axis=1)
out['diff'] = out['PDF'] - out['SmoothedPDF']
return out,tmp_big
if __name__ == "__main__":
out,pchipout = main()
plt.figure()
plt.semilogy(out.index, out.PDF, 'r-', pchipout.index, pchipout.abs().values, 'bo--',
markerfacecolor='none')
plt.legend(['PDF', 'SmoothedPDF (piecewise cubic Hermit spline)'])
plt.figure()
plt.plot(out.index, out.PDF, 'r-', pchipout.index, pchipout.values, 'bo--',
markerfacecolor='none')
plt.legend(['PDF', 'SmoothedPDF (piecewise cubic Hermit spline)'])
plt.show()
Подтверждение стабильности:
>>> out['diff'].any()
False
Логарифмический график для шкалы ошибок:
Как видите, соглашение практически идеально.Обратите внимание на явное отсутствие значений, где ваши исходные данные равны 0. Какой шум остается ниже значений 1e-100 в диапазоне малых х, с которым вам придется учиться жить.
Фактический график интерполированного результата:
