Подгонка данных в пределах двух функциональных границ

Question

У меня есть набор данных, который, согласно теории, ограничен двумя границами.

Верхняя граница: f (x) = 1 / x (x <0.5) </p>

Нижняя граница: f (x) = 1 + 1 / (2x) (x <0,5) </p>

Данные, которые я получил, близки к одной из этих двух границ. Я пытаюсь найти функцию для описания этих данных. Но если я использую нормальный метод подбора, то подгоночная кривая может находиться за пределами этих двух границ. Как я могу заставить свою подогнанную кривую между этими двумя границами, используя scipy.curve_fit? Я пытаюсь использовать аппроксимацию Паде для подгонки, код, который я использую, следующий:

def pfuncp_3_1(x, a0, a1, a2, a3, b1):
    p1 = ((a0+a1*x+a2*x**2+a3*x**3)*(1+b1*x)**(-1)-1-1/(2*x)<0)*2.0
    p2 = ((a0+a1*x+a2*x**2+a3*x**3)*(1+b1*x)**(-1)-1/x>0)*2.0
    return (a0+a1*x+a2*x**2+a3*x**3)*(1+b1*x)**(-1) + p1 + p2


def ffuncp_3_1(x, a0, a1, a2, a3, b1):
    return (a0+a1*x+a2*x**2+a3*x**3)*(1+b1*x)**(-1)


def data_fitter(pfunc, ffunc, fit_x, fit_y, new_x):
    popt, pcov = opt.curve_fit(pfunc, fit_x, fit_ny)
    perr = np.sqrt(np.diag(pcov))
    interp_y = np.zeros(len(interp_x))
    for k in range(len(interp_x)):
        interp_y[k] = ffunc(interp_x[k], *popt)
    return interp_y


if __name__ == "__main__":
    results = data_fitter(pfunc, ffunc, original_x, original_y, interpx)

Answer 1

Во-первых, нужно найти функцию, которую можно отправить на scipy.curve_fit, с некоторыми параметрами, которые могли бы описать нужные нам функции, а именно 1/x и 1+1/(2x): h(x)=a+1/(b*x).

Это даст верхнюю и нижнюю границы для параметров следующим образом:

1/x для a=0, b=1

и

1+1/(2x) для a=1, b=2

Это означает, что a в [0,1] и b в [1,2].

Теперь я смогу сказать scipy.curve_fit, что он должен выбирать параметры a, b и c, чтобы оставаться в пределах:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Generate data
a_true,b_true=.45,1.35
h=lambda x,a,b: a+1/(b*x)
N=100
xmin=0.1
xmax=0.499
x=np.linspace(xmin,xmax,N)
y=h(x,a_true,b_true)+np.random.normal(0,0.05,N)
plt.plot(x,1/x,label="1/x")
plt.plot(x,1+1/(2*x), label="1+1/(2x)")
plt.plot(x,y)
plt.legend()
plt.show()
# Do curve fit
popt, pcov =curve_fit(h,x,y,bounds=[(0,1),(1,2)])
print(popt)
plt.plot(x,1/x,label="1/x")
plt.plot(x,1+1/(2*x), label="1+1/(2x)")
plt.plot(x,h(x,*popt))
plt.show()