[FIXED] Curve fitting with known coefficients in Python

Issue

I tried using Numpy, Scipy and Scikitlearn, but couldn't find what I need in any of them, basically I need to fit a curve to a dataset, but restricting some of the coefficients to known values, I found how to do it in MATLAB, using fittype, but couldn't do it in python.

In my case I have a dataset of X and Y and I need to find the best fitting curve, I know it's a polynomial of second degree (ax^2 + bx + c) and I know it's values of b and c, so I just needed it to find the value of a.

The solution I found in MATLAB was https://www.mathworks.com/matlabcentral/answers/216688-constraining-polyfit-with-known-coefficients which is the same problem as mine, but with the difference that their polynomial was of degree 5th, how could I do something similar in python?

To add some info: I need to fit a curve to a dataset, so things like scipy.optimize.curve_fit that expects a function won't work (at least as far as I tried).

Solution

The tools you have available usually expect functions only inputting their parameters (a being the only unknown in your case), or inputting their parameters and some data (a, x, and y in your case).

Scipy's curve-fit handles that use-case just fine, so long as we hand it a function that it understands. It expects x first and all your parameters as the remaining arguments:

from scipy.optimize import curve_fit
import numpy as np

b = 0
c = 0
def f(x, a):
    return c+x*(b+x*a)

x = np.linspace(-5, 5)
y = x**2

# params == [1.]
params, _ = curve_fit(f, x, y)

Alternatively you can reach for your favorite minimization routine. The difference here is that you manually construct the error function so that it only inputs the parameters you care about, and then you don't need to provide that data to scipy.

from scipy.optimize import minimize
import numpy as np

b = 0
c = 0

x = np.linspace(-5, 5)
y = x**2

def error(a):
    prediction = c+x*(b+x*a)
    return np.linalg.norm(prediction-y)/len(prediction)**.5

result = minimize(error, np.array([42.]))
assert result.success

# params == [1.]
params = result.x

I don't think scipy has a partially applied polynomial fit function built-in, but you could use either of the above ideas to easily build one yourself if you do that kind of thing a lot.

from scipy.optimize import curve_fit
import numpy as np

def polyfit(coefs, x, y):
    # build a mapping from null coefficient locations to locations in the function
    # coefficients we're passing to curve_fit
    # 
    # idx[j]==i means that unknown_coefs[i] belongs in coefs[j]
    _tmp = [i for i,c in enumerate(coefs) if c is None]
    idx = {j:i for i,j in enumerate(_tmp)}

    def f(x, *unknown_coefs):
        # create the entire polynomial's coefficients by filling in the unknown
        # values in the right places, using the aforementioned mapping
        p = [(unknown_coefs[idx[i]] if c is None else c) for i,c in enumerate(coefs)]
        return np.polyval(p, x)
    
    # we're passing an initial value just so that scipy knows how many parameters
    # to use
    params, _ = curve_fit(f, x, y, np.zeros((sum(c is None for c in coefs),)))

    # return all the polynomial's coefficients, not just the few we just discovered
    return np.array([(params[idx[i]] if c is None else c) for i,c in enumerate(coefs)])

x = np.linspace(-5, 5)
y = x**2

# (unknown)x^2 + 1x + 0
# params == [1, 0, 0.]
params = fit([None, 0, 0], x, y)

Similar features exist in nearly every mainstream scientific library; you just might need to reshape your problem a bit to frame it in terms of the available primitives.

Answered By - Hans Musgrave