时间戳

首页 > 解决方案 > 使用边界时如何优化我的鼻窦拟合？

python - 使用边界时如何优化我的鼻窦拟合？

问题描述

我正在尝试对曲线的给定部分进行正弦拟合。有两个限制。首先，我拟合的正弦曲线的偏移量应为 0，其次，我的拟合正弦曲线的幅度应与原始数据的最小值相同。

当我在下面使用我的代码时，配件看起来像我添加的图片（1）。在我看来，窦功能的周期应该更高。拟合曲线至少与我的原始数据匹配，拟合曲线不够宽。当我不使用 c 和 A 的边界时，我的拟合看起来不错（2）。我究竟做错了什么？有没有办法修改拟合，以便在使用 A 和 c 的边界时，正弦曲线更适合？

拟合无界

有界拟合

编辑：我提到的事情是，拟合非常依赖于幅度的起始值（ff_guess）。当我手动将其更改为 X（例如 10 或 30）时，拟合的正弦曲线始终显示接近 X（10.3 或 32.5）的幅度。是否有任何我尚未考虑的设置阻止优化器改变幅度值？

import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from matplotlib.pyplot import figure, rcParams
import numpy as np
from scipy.optimize import curve_fit


#Time
t = [313.544, 313.545, 313.546, 313.547, 313.548, 313.549, 313.55,  313.551, 313.552, 313.553, 313.554, 313.555, 313.556, 313.557, 313.558, 313.559, 313.56,  313.561, 313.562, 313.563, 313.564, 313.565, 313.566, 313.567,]

#y-Values
s = [0.911188, -0.43135, -1.80997, -3.27816, -4.85784, -6.59428, -8.2214, -9.53617, -10.6892, -11.6003, -12.0844, -12.0524, -11.9749, -11.4891, -10.6131, -9.49873, -8.1154, -6.41442, -5.09357, -3.99165, -2.72991, -1.71446, -0.56306, 0.440741]
   
    
#fourier frequency
ff = np.fft.fftfreq(len(t), (t[1]-t[0]))                               
#fourier amplitude
fa = abs(np.fft.fft(s, len(t)))                                            
#Position of maximum Amplitude 
pos_amax = np.argmax(fa[1:])+1                                        
#Frequency at maximum Amplitude (w/2pi)
ff_max = abs(ff[pos_amax])                                             
ff_guess = ff_max   
T_guess = 1000/ff_max 
#A_guess = np.std(s) *2. **0.5                                          
A_guess = min(s)
#c_guess = np.mean(s)
c_guess = 0
#First Guess for all paramters                           
f_guess = np.array([A_guess, 2*np.pi*ff_guess, 0., c_guess])  
#Sinus_Curve
def sin_func(t, A, w, phi, c):
   return A * np.sin(w*t + phi) + c    
#Defining Bounds for A and c
c_bound = 0.1
A_bound = min(s)
#Bounds Array for curve_fit
param_bounds=([1.01*A_bound, -np.inf, -np.inf, -1*c_bound],[0.99*A_bound, np.inf, np.inf, c_bound])
popt, pcov = curve_fit(sin_func, t, s, p0=f_guess, bounds=param_bounds, maxfev=10000000)
#popt, pcov = curve_fit(sin_func, t, s, p0=f_guess, maxfev=10000000)
#
A, w, phi, c = popt 
f = w/(2.*np.pi)
T = 1000/f

t = np.array(t)  
s = np.array(s)

plt.figure(1)
#Generate Sinus Function
s_fit = A * np.sin(w*t + phi) + c
#Plotting
rcParams['figure.figsize'] =10, 5
fig, ax = plt.subplots()  
plt.plot(t, s, "b", label="Original")     
plt.plot(t, s_fit, 'k--', label="Fitting")  
ytitle='ytitle'
xtitle='xtitle'
ax.set(xlabel=xtitle, ylabel=ytitle)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.2f'))               
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))  
ax.grid()
#Sidetext
ausgabe = ("Sinus-Fit \nAmplitude = {:.2f} m/s^2 \nPeriode = {:.2f} ms \nOffset = {:.2f} m/s^2".format(A, abs(T), c))
plt.text(0.795, 0.7, ausgabe, family="sans-serif", fontsize=10, ha='left', va='top', transform=fig.transFigure)  
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.85, box.height]) 
plt.show()

标签： pythonnumpyfftcurve-fitting

解决方案

如果幅度是固定的并且偏移量应该为零，那么为什么要首先安装它们。此外，不需要fft来估计参数，因为有一个简单的线性方法。

看起来像这样：


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


#Time
t = np.array([
    313.544, 313.545, 313.546, 313.547, 313.548, 313.549, 313.55,
    313.551, 313.552, 313.553, 313.554, 313.555, 313.556, 313.557,
    313.558, 313.559, 313.56,  313.561, 313.562, 313.563, 313.564,
    313.565, 313.566, 313.567
])
# ~t -= 313.544

#y-Values
s = np.array([
    0.911188, -0.43135, -1.80997, -3.27816, -4.85784, -6.59428,
    -8.2214, -9.53617, -10.6892, -11.6003, -12.0844, -12.0524,
    -11.9749, -11.4891, -10.6131, -9.49873, -8.1154, -6.41442,
    -5.09357, -3.99165, -2.72991, -1.71446, -0.56306, 0.440741
])

###fixing the amplitude the the maximum value of the data
amp = np.max( np.abs( s ) )

### sine function with the fixed amplitude from above and no offset
def sine(t, w, f ):
    return amp * np.sin( w * t + f )


### getting the nonlinear parameters by using the fact that
### int int y = -y/w**2 + const1 * t + const2
### so the integro-equation as w hidden in a linear factor
### and we can use linear optimization to get it

Sy = cumtrapz( s, x=t, initial = 0 )    ### single integration
SSy = cumtrapz( Sy, x=t, initial = 0 )  ### double integratiom
VMXT = np.array( [ s, t, np.ones( len( t ) ) ] ) ### matrix describing the linear relationship
VMX = np.transpose( VMXT )
A = np.dot( VMXT, VMX )
SV = np.dot( VMXT, SSy )
AI = np.linalg.inv( A )
alpha = np.dot( AI , SV ) ### solution
wstart = np.sqrt( -1 / alpha[0] ) ### as mentioned above, the linear factor ist -1/w**2...so inverse function

### now getting phi by actiully fitting A sin + B cos
### and using tan phi = B/A
VMXT = np.array( [ np.sin( wstart * t ), np.cos( wstart * t ) ] )
VMX = np.transpose( VMXT )
A = np.dot( VMXT, VMX )
SV = np.dot( VMXT, s )
AI = np.linalg.inv( A )
alpha = np.dot( AI , SV ) ### solution
phistart = np.arctan2( alpha[1], alpha[0] )


op, cv = curve_fit( sine, t, s, p0=( wstart, phistart ) )
print( op )
yth = sine( t, *op )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( t, s )
ax.plot( t, yth )
plt.show()

工作得很好。

推荐阅读