Estimating parameter values using optimize.curve.fit

Question

I am attempting to estimate the parameters of the non-linear equation:

y(x1, x2) = x1 / A + Bx1 + Cx2 

using the method outlined in the answer to this question, but can find no documentation on how to pass multiple independent variables to the curve_fit function appropriately.

Specifically, I am attempting to estimate plant biomass (y) on the basis of the density of the plant (x1), and the density of a competitor (x2). I have three exponential equations (of the form y = a[1-exp(-b*x1)]) for the the relationship between plant density and plant biomass, with different parameter values for three competitor densities:

For x2 == 146: y = 1697 * [1 - exp(-0.010 * x1)]

For x2 == 112: y = 1994 * [1 - exp(-0.023 * x1)]

For x2 == 127: y = 1022 * [1 - exp(-0.008 * x1)]

I would therefore like to write code along the lines of:

def model_func(self, x_vals, A, B, C):
    return x_vals[0] / (A + B * x_vals[0] + C * x_vals[1])

def fit_nonlinear(self, d, y):
    opt_parms, parm_cov = sp.optimize.curve_fit(self.model_func, [x1, x2], y, p0 = (0.2, 0.004, 0.007), maxfev=10000)
    A, B, C = opt_parms
    return A, B, C

However I have not found any documentation on how to format the argument y (passed to fit_nonlinear) to capture to two-dimensional nature of the x_vals (the documentation for curve_fit states y should be an N-length sequence). Is what I am attempting possible with curve_fit?

Answer 1

Based on your comment above, you want to think of using the flat versions of the matrices. If you take the same element from both the X1 and X2 matrices, that pair of values has a corresponding y-value. Here's a minimal example

import numpy as np
import scipy as sp
import scipy.optimize

x1 = np.linspace(-1, 1)
x2 = np.linspace(-1, 1)
X1, X2 = np.meshgrid(x1, x2)

def func(X, A, B, C):
    X1, X2 = X
    return X1 / (A + B * X1 + C * X2)

# generate some noisy points corresponding to a set of parameter values
p_ref = [0.15, 0.001, 0.05]
Yref = func([X1, X2], *p_ref)
std = Yref.std()
Y = Yref + np.random.normal(scale=0.1 * std, size=Yref.shape)

# fit a curve to the noisy points
p0 = (0.2, 0.004, 0.007)
p, cov = sp.optimize.curve_fit(func, [X1.flat, X2.flat], Y.flat, p0=p0)

# if the parameters from the fit are close to the ones used 
# to generate the noisy points, we succeeded
print p_ref
print p