0

我正在尝试拟合由反正切函数组成的双断轮廓函数。我的代码似乎不起作用：

XX=np.linspace(7.5,9.5,16)
YY=np.asarray([7,7,7,7.1,7.3,7.5,8.4,9,9.3,9.6,10.3,10.2,10.4,10.5,10.5,10.5])

def func_arc(x,a1,a2,a3,b1,b2,b3,H1,H2):
    beta=0.001
    w1=np.zeros(len(x))
    w2=np.zeros(len(x))
    for i in np.arange(0,len(x)):
        w1[i]=(((math.pi/2)+atan((x[i]-H1)/beta))/math.pi)
        w2[i]=(((math.pi/2)+atan((x[i]-H2)/beta))/math.pi)
    y=(a1*x[i]+b1)*(1-w1[i])+(a2*x[i]+b2)*w1[i]*(1-w2[i])+(a3*x+b3)*w2[i]
    return(y)

其中a和b项是线性回归的斜率和零点值。这些w术语用于切换域。

我考虑了以下连续性限制（H1 y H2）并限制参数：

mask=(XX<=8.2)
mask2=(XX>8.2) & (XX<9)
mask3=(XX>=9)

l1=np.polyfit(XX[mask], YY[mask], 1)
l2=np.polyfit(XX[mask2], YY[mask2], 1)
l3=np.polyfit(XX[mask3], YY[mask3], 1)
H1=(l2[1]-l1[1])/(l1[0]-l2[0])
H2=(l3[1]-l2[1])/(l2[0]-l3[0])

p0=[l1[0],l2[0],l3[0],l1[1],l2[1],l3[1],H1,H2]

popt_arc1, pcov_arc1 =curve_fit(func_arc, XX, YY,p0)

我得到一条线而不是断线（S 形）。
我得到的：

图像显示当前结果

4

1 回答 1

0

这是我的版本。由于线性函数应该连续连接，因此参数实际上更少。因此，偏移量b2和b3没有拟合，而是重新获取线性函数以在过渡处相遇的结果。此外，有人可能会争辩说，数据并不能证明开始或结束时的斜率是合理的。这可以通过减少卡方或其他统计方法来证明/检查。

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

XX=np.linspace( 7.5, 9.5, 16 )
YY=np.asarray( [
    7, 7, 7, 7.1, 7.3, 7.5, 8.4, 9, 9.3, 9.6,
    10.3, 10.2, 10.4, 10.5, 10.5, 10.5
])


yy0 = np.median( YY )
xx0 = np.median( XX )
h0 = 0.8 * min( XX ) + 0.2 * max( XX )
h1 = 0.8 * max( XX ) + 0.2 * min( XX )

def transition( x, x0, s ):
    return ( 0.5 * np.pi + np.arctan( ( x - x0 ) * s ) ) / np.pi

def box( x, x0, x1, s ):
    return transition( x, x0, s ) * transition( -x, -x1, s )

def my_piecewise( x, a1, a2, a3, b1, x0,  x1 ):
    S = 100
    b2 = ( a1 - a2 ) * x0 + b1 ### due to continuity
    b3 = ( a2 - a3 ) * x1 + b2 ### due to continuity
    out = transition( -x , -x0, S ) * ( a1 * x + b1 )
    out += box( x , x0, x1 , S ) * ( a2 * x + b2 )
    out += transition( x , x1, S ) * ( a3 * x + b3 )
    return out

def parameter_reduced( x, a2, b1, x0,  x1 ):
    return my_piecewise(x, 0, a2, 0, b1, x0,  x1 )

def alternative( x, x0, a, s, p, y0 ):
    out = np.arctan( np.abs( ( s * ( x - x0 ) ) )**p )**( 1.0 / p )
    out *= a * (x - x0 ) / np.abs( x - x0 )
    out += y0
    return out

xl=np.linspace( 7.2, 10, 150 )

sol, _ = curve_fit(
    my_piecewise, XX, YY, p0=[ 0, 1, 0, min( YY ), h0,  h1 ] 
)
fl = np.fromiter( ( my_piecewise(x, *sol ) for x in xl ), np.float )

rcp =  np.fromiter( ( y - my_piecewise(x, *sol ) for x,y in zip( XX, YY ) ), np.float )
rcp = sum( rcp**2 )

solr, _ = curve_fit(
    parameter_reduced, XX, YY, p0=[ 1, min(YY), h0,  h1 ]
)
rl = np.fromiter( ( parameter_reduced( x, *solr ) for x in xl ), np.float )

rcpr =  np.fromiter( ( y - parameter_reduced(x, *solr ) for x, y in zip( XX, YY ) ), np.float )
rcpr = sum( rcpr**2 )

guessa = [ xx0, max(YY)-min(YY), 1, 1, yy0 ]
sola, _ = curve_fit(  alternative, XX, YY, p0=guessa)
al = np.fromiter( ( alternative( x, *sola ) for x in xl ), np.float )

rca =  np.fromiter( ( y - alternative(x, *sola ) for x, y in zip( XX, YY ) ), np.float )
rca = sum( rca**2 )

print rcp, rcp / ( len(XX) - 6 )
print rcpr, rcpr / ( len(XX) - 4 )
print rca, rca / ( len(XX) - 5 )

fig = plt.figure( figsize=( 12, 8 ) )
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( XX, YY , ls='', marker='o')
ax.plot( xl, fl )
ax.plot( xl, rl )
ax.plot( xl, al )
plt.show()

结果对我来说看起来不错。

适合

于 2020-09-17T09:28:49.497 回答