Fit data within two function bounds

I have a set of data which, according to a theory, is bounded within two bounds.

The upper bound is : f(x)=1/x (x<0.5)

The lower bound is : f(x)=1+1/(2x)(x<0.5)

The data I got is close to one of these two bounds. I’m trying to find a function to describe these data. But if I use normal fitting method, the fitted curve can be outside of these two bounds. How can I force my fitted curve between these two bounds by using scipy.curve_fit? I'm trying to use Pade approximant to do the fitting, the code I'm using is following:

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)


Show us the code you have so far.
– John Zwinck
Jun 30 at 1:06

The code is added. I use the penalty function to do the fitting@JohnZwinck
– Sizhe
Jun 30 at 1:23

Your question is pretty unclear: Is the data two dimensional, or is it from a one dimensional function? And what do you mean by implementing the bounds? Do you want to obtain them or do you want to find the function?
– Andreas Storvik Strauman
Jun 30 at 9:21

Thanks for your comments @AndreasStorvikStrauman the problem is rephrased
– Sizhe
Jun 30 at 9:34

Aha! I think I get it now.
– Andreas Storvik Strauman
Jun 30 at 9:35

1 Answer
1

First, one needs to find a function that can be sent to scipy.curve_fit, with some parameters, that would be able to describe the functions we want, namely 1/x and 1+1/(2x):
h(x)=a+1/(b*x).

scipy.curve_fit

1/x

1+1/(2x)

h(x)=a+1/(b*x)

This will then give the upper and lower bounds for the parameters as follows:

1/x for a=0,b=1

1/x

a=0

b=1

and

1+1/(2x) for a=1, b=2

1+1/(2x)

a=1

b=2

That means that a is in [0,1] and b is in [1,2].

a

[0,1]

b

[1,2]

Now I'll be able to tell scipy.curve_fit that it has to choose the parameters a,b and c to stay within the bounds:

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()


Thanks, excellent explanation!
– Sizhe
Jun 30 at 13:23

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.