LMFIT：使用多项式模型时约束输出

问题描述

我正在使用 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,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]

解决方法

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

        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')

,

一种选择可能是使用样条曲线。 一个快速而肮脏的方法，只是为了展示这个想法，可能看起来像这样：

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,.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 ) )


## 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 样条来建立不动点为固定值，即不拟合其参数。