Constraining the model value at a certain point by using an expression

[fmh...@gmail.com]

Hi there,

I am trying to find the best set of parameters for each of several decay experiments. Most importantly, I need those parameters to produce a model, that attains a prescribed value y_m at a given location x_m. Only as a second-most important requirement, I want a good fit.

The following image illustrates my current status: The fit is good, except at the maximum, where I need the fit to attain the same value as the data,

i.e. (x_m, y_m) = (1, 1).

I tried to achieve that by using an expression as in the following example:

def weibull(x, k, lam, amp, stat):
    return (amp - stat) * k / lam * (x / lam) ** (k - 1) * np.exp(-(x / lam) ** k) + stat

mod = lmfit.Model(weibull)
pars = lmfit.Parameters()
pars.add('k', value=k0, vary=True, min=0.0)
pars.add('lam', value=1, vary=True, min=0.0)
pars.add('amp', value=1, vary=True)
pars.add('stat', value=0.08, vary=False)

pars._asteval.symtable['fct'] = weibull
pars.add('constraint', expr='fct(1.0, k, lam, amp, stat)', value=1.0, vary=False)

I don't get any error messages, but the returned value of the parameter 'constraint' differs largely from the prescribed value 1. In other words, the resulting fit does not change at all by including the constraint. I am guessing, that this usage is not supported. Is there an alternative approach to what I want to achieve within lmfit? If that requirement goes beyond what is possible with lmfit, do you have an idea about other available algorithms/approaches that could help me?

Any help is very appreciated!

Thank's so far and best wishes,

Florian

[Matt Newville]

It appears that you want to have your model (say, 'y') take some specific value for some specific value of the independent variables (say, 'x').   

That is not exactly what constraint expressions do.  They constrain the parameters that are used to calculate the model.  The reason the fitting does not change when adding your constraint expression is because you do not use the `constraint` parameter in your model function. 

What you might intend to do is

    def weibull(x, k, lam, amp, stat):

        model = (amp - stat) * k / lam * (x / lam) ** (k - 1) * np.exp(-(x / lam) ** k) + stat

        model[np.where(x==1.0)]  = 1.0

        return model

But, I'm not sure that really is a great idea.   Probably a better approach is really to use a constraint on one of the parameters (`stat`, `amp`, `lam`,  `k`) so that the value of the function is 1.0 at x=1. I think you should be able to work out an expression for `amp` in terms of the parameters `k`, `stat`, and `lam` that guarantees model[x=1] = 1.

--Matt

[fmh...@gmail.com]

Hi Matt,

thanks for your quick answer!

Yea, that makes sense. I'll try to come up with an expression for 'amp'..

Best wishes,

Florian