0

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?

Community
  • 1
  • 1
shbrainard
  • 377
  • 2
  • 9
  • I don't quite understand your question. Is `y(x1, x2)` a scalar function (i.e., given two numbers x1 and x2 it returns a single number? If so, then the argument `y` passed to `fit_nonlinear` shouldn't have to be two-dimensional. If `y(x1, x2)` is supposed to be a vector function (it returns two numbers for each x1, x2 pair), then you can split it up into two scalar functions (`y1(x1,x2)` and `y2(x1,x2)`) and fit them each individually. – hunse Dec 02 '13 at 21:36
  • The former: y(x1,x2) returns a single number given a value for x1 and x2. However, I am not sure how to structure the y-values as an N-length sequence so that their relationship to the 2-dimensional matrix of x-values is clear. – shbrainard Dec 03 '13 at 11:00

1 Answers1

0

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
hunse
  • 3,175
  • 20
  • 25