来自已知系数的 Scipy BSpline 对象

问题描述

给定节点 x_i 处的观测值 y_i，在 x_i 处找到基函数 B_j(x) 及其系数 c_j 的最小二乘优化问题由矩阵方程给出

c = (B'*B+ lam * D'*D)^(-1)*B'y

B_ij 是在 x_i 处评估的 BSplines 矩阵，D_ij 是惩罚导数粗糙度的差分矩阵。有关详细信息，请参阅 here 或此 crash course。

如果 lam=0，则可以求解该方程组，使得 B*c=y。事实上，我对 BSplines 的实现正是给出了这个结果。但是，如果我现在使用 BSplines 的 Scipy 实现，则

spl=BSpline(x,c,3)

边缘失败。这是怎么回事？

def BasisFunction(t,i):
   if np.abs(t-i)<1:
       return (4/6)-np.abs(t-i)**2+0.5*np.abs(t-i)**3
   elif np.abs(t-i)<2:
       return (1/6)*(2-np.abs(t-i))**3
   else:
       return 0
   
   
def GetSpline(x,y,lam):
   if len(x)>8:
       #tic = time.time()
       x=np.array(x)
       y=np.array(y)
       #print(x)
       dx=min(np.diff(x))
       x_EquallySpaced=np.arange(x[0],x[-1]+dx,dx) #Prepare for derivative calculation
       y_EquallySpaced=sc.interpolate.splev(x_EquallySpaced,sc.interpolate.splrep(x,y)) #Piecewise Linear Interpolation for derivative
       #print(y_EquallySpaced)
       
       B = np.transpose([[BasisFunction(x_EquallySpaced[j],x_EquallySpaced[i]) for i in range(0,len(x_EquallySpaced))] for j in range(0,len(x_EquallySpaced))])
       print(B)
       D=(1/dx**3)*sc.sparse.diags([-0.5,1,-1,0.5],[-2,2],shape=np.shape(B)).toarray() #Third Derivative Finite Central Difference matrix
       D[0:2,:]=0
       D[-2:,]=0
       

       #print(np.matmul(np.transpose(D),D))
       coefs =sc.linalg.pinv(np.transpose(B)@B+lam*np.transpose(D)@D)@np.transpose(B)@y_EquallySpaced

   return [B,coefs,BSpline(x_EquallySpaced,3,extrapolate=True)]

plt.figure()
#plt.subplot(2,1)
plt.plot(Data['time'].loc[0],res[2](Data['time'].loc[0]),label='Scipy BSpline')
plt.plot(Data['time'].loc[0],res[0]@res[1],label='Direct')
plt.plot(Data['time'].loc[0],Data['x'].loc[0],'x',label='Data')
plt.legend()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

bspline interpolation scipy