LMFIT：使用多项式模型时限制输出答案

【问题标题】：LMFIT: Constraining the output when using the polynomial modelLMFIT：使用多项式模型时限制输出
【发布时间】：2021-05-31 23:19:35
【问题描述】：

我正在使用 LMFIT 将分段多项式拟合到正弦波的第一象限。我希望能够在多项式输出上添加一个约束——而不是对其参数。例如，我想确保输出是 >= 0 和

我知道使用 np.polyfit 可能会更好，但最终我想添加更多非线性约束，并且 LMFIT 框架更灵活。

import numpy as np
from lmfit.models import LinearModel

#split sine wave in 4 segments with 1024 points
nseg = 4
frac = 2**10
npoints = nseg*frac
xfrac = np.linspace(0, 1, num=frac, endpoint=False)
x = np.linspace(0, 1, num=npoints, endpoint=False)
y = np.sin(x*np.pi/2)
yseg = np.reshape(y, (nseg, frac))


mod = LinearModel()

coeff = []
bestfit = []
for i in range(nseg):
    pars = mod.guess(yseg[i], x=xfrac)
    out = mod.fit(yseg[i], pars, x=xfrac)
    coeff.append([out.best_values['slope'], out.best_values['intercept']])
    bestfit.append(out.best_fit)
bestfit = np.reshape(bestfit, (1, npoints))[0]

【问题讨论】：

我只能看到如何在模型参数上添加约束，而不是它的输出（这里是每个段中 y=mx+b 的结果）
我是否正确地说您的分段函数不是连续的？关于输出参数的限制，我猜你必须手动计算输入参数的相应限制。
你说得对，它不是连续的。强制连续性似乎不是一个通用的约束，所以我认为需要完成与下面发布的相同的想法。
你可以尝试线性 B 样条的意义

标签： python constraints curve-fitting lmfit piecewise

【解决方案1】：

事实证明，这是通过对参数本身添加约束来完成的，这些约束变成了模型输出的正确约束。使用自定义模型进行线性插值可以如下完成：

        def func(x, c0, c1):
            return c0 + c1*x
        pmodel = Model(func)
        params = Parameters()
        params.add('c0')
        params.add('clip', value=0, max=1.0, vary=True)
        params.add('c1', expr='clip-c0')

【讨论】：

【解决方案2】：

一个选项可能是使用样条曲线。一个快速而肮脏的方法，只是为了呈现这个想法，可能是这样的：


import matplotlib.pyplot as plt

import numpy as np

## quich and dirty spline function
def l_spline(x, abc ):
    if isinstance( x, ( list, tuple, np.ndarray ) ):
        out = [ l_spline( elem, abc ) for elem in x]
    else:
        a, b, c = abc
        if x < a:
            f = lambda t: 0
        elif x < b:
            f = lambda t: ( t - a ) / ( b - a )
        elif x < c:
            f = lambda t: -( t - c ) / (c - b )
        else:
            f = lambda t: 0
        out = f(x)
    return out

### test data
xl = np.linspace( 0, 4, 150 )
sl = np.fromiter( ( np.sin( elem ) for elem in xl ), np.float )

### test splines with manual double knots on first and last
yl = dict()
yl[0] = l_spline( xl, ( 0, 0, .4 ) )
for i in range(1, 10 ):
    yl[i] = l_spline( xl, ( (i - 1 ) * 0.4 , i * 0.4, (i + 1 ) * 0.4 ) )
yl[10] = l_spline( xl, ( 3.6, 4, 4 ) )


## This is the most simple linear least square for the coefficients

AT = list()
for i in range( 11 ):
    AT.append( yl[i] )
AT = np.array( AT )
A = np.transpose( AT )

U = np.dot( AT, A )
UI = np.linalg.inv( U )
K = np.dot( UI, AT )
v = np.dot( K, sl )


## adding up the weigthed sum
out = np.zeros( len( sl ) )
for a, l in zip( v, AT ):
    out += a * l

### plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( xl, sl, ls=':' )
for i in range( 11 ):
    ax.plot( xl, yl[i] )
ax.plot( xl, out, color='k')
plt.show()

看起来像这样：

而不是简单的线性优化，可以使用更复杂的功能来确保没有参数大于1.这会自动确保功能不超过1.可以通过设置根据B样条来建立一个固定点到固定值，即不符合其参数。

【讨论】：