4

Is there a way to construct a an lmfit Model based on a function with an arbitrary number of dependent variables? For example:

from lmfit import Model

def my_poly(x, *params):
  func = 0
  for i in range(len(params)):
    func+= params[i]*z**i
  return func

#note: below does not work
my_model = Model(my_poly, independent_vars = ['x'], param_names = ['A','B','C'])

Something similar to the above would be wonderful if I am interested in a polynomial series and want to test the performance as the series grows or shrinks.

Brian Pollack
  • 198
  • 3
  • 12

1 Answers1

4

Since Model() uses function argument names to build parameter names, using *params won't work easily (how would one know to call them A, B, C, and not coeff0, coeff1, coeff2, or something else?).

I don't know that a truly arbitrary number could be supported, but it should be possible to do a very large number. The Polynomial Model (see http://lmfit.github.io/lmfit-py/builtin_models.html#polynomialmodel and https://github.com/lmfit/lmfit-py/blob/master/lmfit/models.py#L126 for implementation) supports up to 7 coefficients. It should be no problem to extend that to a much larger number. It might easily lead to computational problems, but I think that is what you are expecting to explore.

If you're willing to make a small change, it is be possible to do something like you're looking for. This uses keyword arguments instead of positional arguments, and relies on parameter name order (that is with sort) to indicate which coefficient goes with what exponent, rather than order of the positional arguments. This might be close to what you're looking for:

import numpy as np

from lmfit import Model, Parameters

def my_poly(x, **params):
    val= 0.0
    parnames = sorted(params.keys())
    for i, pname in enumerate(parnames):
        val += params[pname]*x**i
    return val

my_model = Model(my_poly)

# Parameter names and starting values
params = Parameters()
params.add('C00', value=-10)
params.add('C01', value=  5)
params.add('C02', value=  1)
params.add('C03', value=  0)
params.add('C04', value=  0)

x = np.linspace(-20, 20, 101)
y = -30.4 + 7.8*x - 0.5*x*x + 0.03 * x**3 + 0.009*x**4
y = y + np.random.normal(size=len(y), scale=0.2)

out = my_model.fit(y, params, x=x)
print(out.fit_report())

Hope that helps.

M Newville
  • 7,486
  • 2
  • 16
  • 29
  • Excellent! I've explored the built-in polynomial model, but it is difficult to use in my actual analysis. However, you suggestion regard kwargs instead of args is perfect. Thank you. – Brian Pollack Sep 08 '15 at 16:02
  • Whille were are on this topic of arbitrary number of parameters, would it be possible to say fit an arbitrary number of gaussians ? i.e where the number of gaussian models is variable ? of course with limits say from 2 - 10 max ? Thanks in advance. – user1301295 Jun 29 '16 at 22:33