How to put a constraint on the expression in a natuarl logarithm function that I am trying to fit into my data?

[nepomnyi]

Hello dear lmfit community.

I am having problems fitting natural logarithm function into my data - the model keeps getting NaNs and aborting. I think this might be happening due to the model yielding a negative expression within the ln() function. Of course there might be other reasons, but for now, I want to make sure the model never gets a negative expression and see if it work.

I have discovered this page that talks about constraints, but it is not helpful for my case. Thus, I tried to infer from the said page how to impose a constraint on the entire expression within the ln() function. Here is what I came up with (remember, the expression within the ln() function must always be positive):

# First, I create the most general form of the ln() function.
def lnFit(coord, p1, p2, p3, p4):
        return p1 * np.log(p2 * coord + p3) + p4

gmodel = Model(lnFit)

params = Parameters()
params.add('p1', value=4.9, min=0.000000001) # I want the first parameter be positive (for my own reasons)
params.add('p4', value=1.1) # I don't want any constraints for it

# Then I am following the referenced page. I want (p2*coord+p3) be positive always .
# I follow the example "x+y<=10" from the referenced page:
params.add('delta', value=0.000000001, min=0.000000001, vary=True) # I want my expression be positive ==> delta must be positive
params.add('coord', value=5, vary=True) # This is the problem: I do not know how to incorporate independent variable into this procedure - I decided to make it a parameter
params.add('p3', value=1.8, min=0.000000001) # I want p3 be positive (for my own reasons)
params.add('p2', expr='(delta-p3)/coord') # this is the expression: p2 can assume only those values that will make p2 * coord + p3 > delta

# Then I fit the model. But I have a problem now: I added parameters that are not in lnFit() function ('delta' and 'coord'). So, I am specifying 
# the parameters manually instead of doing "params=params" below:
result = gmodel.fit(Y, coord=X, p1=params['p1'], p2=params['p2'], p3=params['p3'], p4=params['p4']) # here Y and X are my data

Having run the code, I get an error message: NameError: name "delta" is not defined. And I don't have any other ideas.

How do I make sure that p2 * coord - p3 always stays positive during the fitting procedure?

Thank you in advance.
Ivan

---

UPDATE 1.
I changed the function as given below:

def lnFit(coord, p1, p2, p3, p4):
        return p1 * np.log(p2 * (coord + p3)) + p4

Then using my data I figured out 1) the largest possible value of coord I can have and 2) therefore, the shift parameter p3. Because my ln() function must be "facing" the left (not the right as it normally does), I could figure out the largest possible value of parameter p2. In other words, I just used min and max arguments in Parameters().add() function - I did not do any < or > type constraints. I proceeded fitting the model, it succeeded, but the fit is just wrong.

Irrespectively of all this, my original question still holds: how does one impose a constraint on an expression that has the independent variable in it?

---

UPDATE 2.
Like I mentioned in UPDATE 1 above, I managed to get some fit, but it was just wrong. I couldn't find a way to make it right and abandoned the logarithmic fit altogether. I ended up fitting a sigmoid function instead of the logarithmic function. I used the logistic function for sigmoid and lmfit fitted it from the first attempt without any parameter adjustment from my side.
This is interesting for the following reason. As you may notice the logistic function is based on the exponent. I had another dataset that looked like an exponential decay function which I fitted in without any problems with the help of lmfit. Exponential decay function is also based on the exponent. Moreover, it is the "inverse" of my troublesome logarithmic function. I.e., the pattern that I see is that lmfit has no problems fitting exponent-based functions, but it starts sucking when it comes to logarithms. I am sure this is because of the constraint that we cannot take a logarithm of a negative number (whereas, we can raise the exponent to whatever power we like). This leads me to a suggestion: lmfit should take care of situations when a function cannot be taken of negative numbers.

[newville]

My main advice would be to prevent negative arguments to np.log() - that would be independent of the values of your coord and the parameter values. Like, you could just use an "abs()" or look for and negative values with np.where((p2*coord + p3)<0) in your lnfit function. That strengthens the function.

If coord is all positive or all negative, and if p3 is positive, setting an upper or lower limit on p2 shouldn't be too hard: just p2 > 0 or p2<0 would be sufficient.

But if coords has both positive and negative values, it is much trickier, as the current value of p3 would need to be known to set a bound on p2. And p3 not being known is sort of one of the points. Since you are asking the question, I sort of think that coords be both positive and negative must be a possibility. I think the general solution is simply to ensure that p2*coord+p3 is positive, which brings me back to the original conclusion.

For an inequality constraint (as shown in the example), you would have to define a parameter delta if you are going to use that in an expression. But, I don't think that is really going to help here if you don't know the sign of coords.

[nepomnyi]

Thank you @newville . I did almost exactly what you are suggesting and described that in the update to my post (may be not as articulated as you did). But that's a pretty convoluted way around and, most certainly, this is not a general solution to the problem. And I couldn't make it work anyway (see the second update to my post).

To summarize: lmfit does not have an option of including independent variable into the expressions, is that correct?

Lastly, you mentioned that I will have "to define a parameter delta" if I want to use it. Could you, please, update the referenced page to show how to do it? Right now, the example there is incomplete and I can only guess how to "define" delta. I made a guess as I showed in my post and it was wrong. The referenced page only shows how to add delta to Parameters(). But because delta is not present in my function that I am trying to fit lnFit(coord, p1, p2, p3, p4), lmfit complains and gives me an error. So, if I could use this delta approach, what would be the complete example for the x+y<=10 shown on the referenced page?

[newville]

@nepomnyi

If you are not putting a check of the values sent to np.log() in your model function, then you are not doing what I suggested. I don't think trying to solve this with constraints is going to be the better solution.

To summarize: lmfit does not have an option of including independent variable into the expressions, is that correct?

Um, kind of "neither". Constraints use mathematical expressions within an asteval interpreter. Independent variables can be any python object, and could, in principle, be added to the asteval interpreter. But, you would have to do that, and then use that.

That is a very convoluted appoach to preventing log() from getting a negative value, when you could just prevent log() from getting a negative value.

Lastly, you mentioned that I will have "to define a parameter delta" if I want to use
it. Could you, please, update the referenced page to show how to do it? Right now,
the example there is incomplete and I can only guess how to "define" delta.

The example shows how to define delta.

I made a > guess as I showed in my post and it was wrong. The referenced page only shows how to
add delta to Parameters(). But because delta is not present in my function that I am
trying to fit lnFit(coord, p1, p2, p3, p4), lmfit complains and gives me an error.

Um, no it doesn't... Maybe you are getting an error from your very weird usage of passing parameter values to Model.fit.

result = gmodel.fit(Y, coord=X, p1=params['p1'], p2=params['p2'], p3=params['p3'], p4=params['p4']) # here Y and X are my data

I don't know what this is supposed to do. But you did not post any error messages from that...

So, if I could use this delta approach, what would be the complete example for the
x+y<=10 shown on the referenced page?

Just to move to a completely different problem domain (because I am not going to do your homework) you can do something like:

import numpy as np
from lmfit import Model, Parameters
from lmfit.lineshapes import gaussian

def two_peaks(x, amp1, cen1, sig1, amp2, cen2, sig2):
    return gaussian(x, amp1, cen1, sig1) + gaussian(x, amp2, cen2, sig2)

pars = Parameters()
pars.add('amp1', value=35.0)
pars.add('amp2', value=20.0)
pars.add('sig1', value=6.0, min=0)
pars.add('sig2', expr='2.0*sig1')
pars.add('cen1', value=50.0, min=0, max=100)
pars.add('peak_sep', value=5.0, min=0)
pars.add('cen2', expr='cen1 + peak_sep')

x = np.linspace(0, 100, 501)
y = gaussian(x, 30, 45, 3.5) + gaussian(x, 15, 57, 7.05)
y += np.random.normal(size=len(x), scale=0.1)

model = Model(two_peaks)
ret = model.fit(y, params=pars, x=x)
print(ret.fit_report())

Note that peak_sep is varied in the fit, and is not used directly by the model function. That is OK.
Lmfit will not complain about this or give an error.

[nepomnyi]

Thank you for your response. And if you could add your last example to the docs, that would be greatly appreciated. For instance, I had no idea it would be "OK" like you wrote if I just added peak_sep to Parameters() and then did params=pars like you showed. I was thinking that, because peak_sep is not one of the arguments of two_peaks, I could not use params=pars and had to specify each parameter manually. This is an example of confusion that could be resolved by adding your full example of two-peak to https://lmfit.github.io/lmfit-py/constraints.html.

[newville]

See https://lmfit.github.io/lmfit-py/examples/example_fit_with_inequality.html#sphx-glr-examples-example-fit-with-inequality-py