Constraints does not work in non-linear regression

[thim peng]

Hi all,

I'm doing non-linear regression using some data, and I want the result curve crosses a certain point, so I add some constraints in parameters setting, like this:

# add constrain condition

def params_constraint(init_params,lb,ub,varies):

    params_second_step = Parameters()

#     params_second_step.add('a', value = init_params['a'], min=lb['a'], max=ub['a'],vary=varies['a'])

    params_second_step.add('b', value = init_params['b'], min=lb['b'], max=ub['b'], vary=varies['b'])

    params_second_step.add('c', value = init_params['c'], min=lb['c'], max=ub['c'], vary=varies['c'])

    params_second_step.add('d', value = init_params['d'], min=lb['d'], max=ub['d'], vary=varies['d'])

    params_second_step.add('theta', value = init_params['theta'], min=lb['theta'], max=ub['theta'], vary=varies['theta'])

    params_second_step.add('tb', value = init_params['tb'], min=lb['tb'], max=ub['tb'], vary=varies['tb'])

    params_second_step.add('a', value = init_params['a'], min=lb['a'], max=ub['a'],vary=varies['a'],expr = '(tb-(b*pow(sin(c*theta),2)+pow(cos(c*theta),2))*d)/pow(theta,2)')

#     params_second_step.add('descr', value = 0, vary=False, expr = ('a*pow(theta,2)+(b*pow(sin(c*theta),2)+pow(cos(c*theta),2))*d-tb'))

    return params_second_step

I tried to set constraint parameter 'a' with regression equation and set a parameter 'descr' as fix value, but in some points, the constraint seems can't work.

Any help is very appreciated!

Best

Peng

[Matt Newville]

What does not work?  

[Matt Newville]

On Thu, Feb 6, 2020 at 7:51 AM thim peng <zsj...@gmail.com> wrote:

the constraint conditions sometimes does not work means the result curve can not cross the certain point as I need the regression equation exactly satisfy at the point, for example the red line below should cross the red point below.

Sorry, I have no idea what you're doing or looking for.   Why do you think one curve cannot cross a certain point?   Should that be obvious to me?   It is not.

Your example code is far too messy to read and has far too much extraneous junk.   Make it simpler and clearer, and tell us what your are expecting. 

--Matt

[thim peng]

Hi Matt,

Sorry for my unclear description and messy code, I have made it simpler and clear.

I want do some non-linear regression on a group data (theta1, yh1) with the function defined as fit_function(x,a,b), it has 2 parameters: a and b. 

def fit_function(x,a,b):

    u = np.radians(x)

    return a*np.power(u,2)+(b*np.power(np.sin(u),2)+np.power(np.cos(u),2))*271.4877

What's more, I want the regression result curve just passing through another certain point (theta2, yh2), that means the result of fit_function(theta2,a,b)-yh2 should exactly equals 0.

So I add some constraints on a by setting expression using fit_function, and also try to add another fixed parameter 'descr' : fit_function(theta2,a,b)-yh2 as 0, 

model = Model(fit_function,indepent_vars=['x'],)

params = Parameters()

# params.add('a', value = -1.1, max=0)

params.add('b', value = 0.5, max=1,vary=True)

params.add('theta', value = np.radians(theta2), vary=False)

params.add('tb', value = yh2, vary=False)

params.add('a', value = -1.1, max=0,expr = '(tb-(b*pow(sin(theta),2)+pow(cos(theta),2))*271.4877)/pow(theta,2)')

# params.add('descr', value = 0, vary=False, expr = '(a*pow(theta,2)+(b*pow(sin(theta),2)+pow(cos(theta),2))*271.4877-tb)')

but both of the finnal result of fit_function(theta2,a,b)-yh2 are still not 0.

result = model.fit(yh,x=theta1,params=params, weights=1/xacc)

print('check result should be 0:',fit_function(theta2,result.params['a'].value,result.params['b'].value)-yh2)

and out print is:

check result should be 0: -5.643670823520608

I can get the expected result 0 by scipy, but I wonder why I cannot get the same one by lmfit 1.0.0, do you have any ideas?

Thanks again for your patience and response.

在 2020年2月7日星期五 UTC+8上午2:35:42，Matt Newville写道：

[Omer Luria]

Hi,

I did not get too much into your code, but why don't you just constraint this within your fit_function? It is almost always better to have most information contained in your model function.

In this case it's very easy, just substitute fit_function(theta2)=yh2 and you can express a with b, theta2 yh2 thus you eliminate one of them.

Best,

Omer

[Matt Newville]

On Fri, Feb 7, 2020 at 1:25 AM thim peng <zsj...@gmail.com> wrote:

Hi Matt,

Sorry for my unclear description and messy code, I have made it simpler and clear.

I want do some non-linear regression on a group data (theta1, yh1) with the function defined as fit_function(x,a,b), it has 2 parameters: a and b. 

def fit_function(x,a,b):

    u = np.radians(x)

    return a*np.power(u,2)+(b*np.power(np.sin(u),2)+np.power(np.cos(u),2))*271.4877

What's more, I want the regression result curve just passing through another certain point (theta2, yh2), that means the result of fit_function(theta2,a,b)-yh2 should exactly equals 0.

So I add some constraints on a by setting expression using fit_function, and add another fixed parameter 'descr' : fit_function(theta2,a,b)-yh2 as 0, 

model = Model(fit_function,indepent_vars=['x'],)

params = Parameters()

# params.add('a', value = -1.1, max=0)

params.add('b', value = 0.5, max=1,vary=True)

params.add('theta', value = np.radians(theta2), vary=False)

params.add('tb', value = yh2, vary=False)

params.add('a', value = -1.1, max=0,expr = '(tb-(b*pow(sin(theta),2)+pow(cos(theta),2))*271.4877)/pow(theta,2)')

# params.add('descr', value = 0, vary=False, expr = '(a*pow(theta,2)+(b*pow(sin(theta),2)+pow(cos(theta),2))*271.4877-tb)')

but the finnal result of fit_function(theta2,a,b)-yh2 is still not 0.

result = model.fit(yh,x=theta1,params=params, weights=1/xacc)

print('check result should be 0:',fit_function(theta2,result.params['a'].value,result.params['b'].value)-yh2)

and out print is:

check result should be 0: -5.643670823520608

I can get the expected result 0 by scipy, but I wonder why I cannot get the same one by lmfit 1.0.0, do you have any ideas?

Remove the bounds on `a` and `b` and it will work.   Some general advice:

  a) ALWAYS look at the fit report.  It will, for example, tell you when parameters are at bounds.

  b)  use bounds when there is a mathematical of a physical reason for a parameter not exceeding some value.

  c) use `x**2` instead of `pow(x, 2)`.  Readability counts.

[thim peng]

Hi Omer,

Thanks for your tips very much, and I get the expected result by express a as other paramters.

Thanks again.

Best

Peng

在 2020年2月7日星期五 UTC+8下午5:26:17，Omer Luria写道：

[thim peng]

Hi Matt,

Thanks for your advices very much, and I can get the same result with tips given by Omer when I remove bounds of a, and your advices show me another point of view to my questions.

Thanks again.

Best

Peng

在 2020年2月7日星期五 UTC+8下午8:55:20，Matt Newville写道：