Multivariate curve fit in python

Question

Can somebody please point me in the right direction...

I need to find the parameters a,b,c,d of two functions:

Y1 = ( (a * X1 + b) * p0 + (c * X2 + d) * p1 ) / (a * X1 + b + c * X2 + d)

Y2 = ( (a * X2 + b) * p2 + (c * X2 + d) * p3 ) / (a * X1 + b + c * X2 + d)

X1, X2 (independent variables) and Y1, Y2 (dependent variables) are observations, i.e. one-dimensional arrays with thousands of entries each.

p0, p1, p2, p3 are known constants (scalars).

I successfully solved the problem with the first function only with a curve-fit (see below), but how do i solve the problem for Y1 and Y2 ?

Thank you.

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


X  = [X1,X2]

def fitFunc(X, a,b,c,d):
    X1, X2 = X
    return ((a * X1 + b) * p0  + (c * X2 + d) * p1) / (a * X1 + b + c * X2 + d)

fitPar, fitCov = curve_fit(fitFunc, X, Y1)

print(fitPar)

Answer 1

One way would be to minimize both your functions together using scipy.optimize.minimze. In the example below, a function residual is passed a, b, c, and d as initial guesses. Using these guesses, Y1 and Y2 are evaluated, then the mean squared error is taken using the data and predicted values of respective functions. The error is returned as the mean error of the two functions. The optimized set of parameters is stored in res as res.x.

import numpy as np
from scipy.optimize import minimize

#p0 = ... known
#p1 = ... known
#p2 = ... known
#p3 = ... known

def Y1(X, a,b,c,d):
    X1, X2 = X
    return ((a * X1 + b) * p0  + (c * X2 + d) * p1) / (a * X1 + b + c * X2 + d)

def Y2(X, a,b,c,d):
    X1, X2 = X
    return ((a * X1 + b) * p2  + (c * X2 + d) * p3) / (a * X1 + b + c * X2 + d)

X1 = np.array([X1]) # your X1 array
X2 = np.array([X2]) # your X2 array
X = np.array([X1, X2])
y1_data = np.array([y1_data]) # your y1 data
y2_data = np.array([y2_data]) # your y2 data

def residual(x):
    a = x[0]
    b = x[1]
    c = x[2]
    d = x[3]
    y1_pred = Y1(X,a,b,c,d)
    y2_pred = Y2(X,a,b,c,d)
    err1 = np.mean((y1_data - y1_pred)**2)
    err2 = np.mean((y2_data - y2_pred)**2)
    error = (err1 + err2) / 2
    return error

x0 = [1, 1, 1, 1] # Initial guess for a, b, c, and d respectively
res = minimize(residual, x0, method="Nelder-Mead")  
print(res.x)