CODEWITHSUNDEEP
pythonnumpyscipycurve-fitting
bob.sacamento
bob.sacamento

Reputation: 6651

two dimensional fit with python

I need to fit a function

z(u,v) = C u v^p

That is, I have a two-dimensional data set, and I have to find two parameters, C and p. Is there something in numpy or scipy that can do this in a straightforward manner? I took a look at scipy.optimize.leastsq, but it's not clear to me how I would use it here.

Upvotes: 5

Views: 5795

Answers (3)

David K.
David K.

Reputation: 477

On the risk of being a bit late to the party, I would like to elaborate on the other answers. Flattening does work however from this answer it is not really clear why and honestly more of a bug than a feature.

The package fitting_toolkit provides a function to do this in a straightforward way building on scipy.optimize. The advantage here is that fitting_toolkit.multivariate_fit() does not care at all how you format the data as it only works within the parameter space.

Full disclosure, I'm the lead contributor on that package. I implemented this because I needed it.

First define your model:

def z(uv, C, p):
   return C * uv[0] * uv[1]**p

Where uv is a tuple containing numpy arrays of your u and your v data. How exactly you package this is up to you, the library does not touch your data. Just make sure your model returns a 1- or 0-dimensional array to sum over. You can then fit your model:

from fitting_toolkit import multivariate_fit as fit

uv = (u, v)
sigma = np.ones_like(z_data) # mock error for equal weights
popt, pcov = fit(z, uv, z_data, sigma=sigma, theta_0=theta_0)

Brief explanation for DIY

Fitting a model means minimizing a loss function. You can do that yourself with scipy.optimize.minimize(Docs).

from scipy.optimize import minimize

# Loss function for unweighted least squares
def loss_function(parameters):
    return np.sum((output - model(input, *parameters))**2)

result =  minimize(loss_function, theta_0, **kwargs)

popt = result.x
pcov = result.hess_inv

Upvotes: 0

DavidW
DavidW

Reputation: 30889

def f(x,u,v,z_data):
  C = x[0]
  p = x[1]
  modelled_z = C*u*v**p
  diffs = modelled_z - z_data
  return diffs.flatten() # it expects a 1D array out. 
       # it doesn't matter that it's conceptually 2D, provided flatten it consistently

result = scipy.optimize.leastsq(f,[1.0,1.0], # initial guess at starting point
                        args = (u,v,z_data) # alternatively you can do this with closure variables in f if you like
                              )
# result is the best fit point

For your specific function you might be able to do it better - for example, for any given value of p there is one best value of C that can be determined by straightforward linear algebra.

Upvotes: 6

Frank M
Frank M

Reputation: 1570

You can transform the problem into a simple linear least squares problem, and then you don't need leastsq() at all. 

z[i] == C * u[i] * v[i]**p

becomes

z[i]/u[i] == C * v[i]**p

And then

log(z[i]/u[i]) == log(C) + p * log(v[i])

Change variables and you can solve as a simple linear problem:

Z[i] == L + p * V[i]

Using numpy and assuming you have the data in arrays z, u and v, this is rendered as:

Z = log(z/u)
V = log(v)
p, L = np.polyfit(V, Z, 1)
C = exp(L)

You probably ought to put a try: and except: around it in case some of the u values are zero or there are negative values.

Upvotes: 1

Related Questions