我的受限（自然）三次样条方程错了吗？

计算科学 插值 r b样条
2021-11-27 08:16:35

我正在尝试将 4 节的受限三次样条（自然三次样条）拟合到玩具数据中，尝试遵循 Hastie、Tibshirani、Friedman 2nd ed。5.2.1 p.144-146，方程 5.4 和 5.5。数据：基本上是一个转置的“S”形。R代码是：

n <© 100
x <- (1:n)/n
true <- ((exp(1.2*x)+1.5*sin(7*x))-1)/3
noise <- rnorm(n, 0, 0.15)
y <- true + noise
plot(x,y)

我将结设置为：{.2, .4, .6, .8} 并使用 R 中的非线性 NLS() 函数进行拟合，但无论如何我都无法获得数据的 S 形尝试。

我的方程式错了吗？或者我的方法完全偏离了基础？有什么建议？

（书摘、我的方程式和下面发布的数据图）

在此处输入图像描述

自然三次样条K结表示为K基函数。可以从三次样条的基础开始，并通过施加边界约束来导出简化的基础。例如，从 5.2 节中描述的截断幂级数基础开始，我们得到

dk=(X−ξk)+3−(X−ξK)+3ξK−ξk.
可以看出这些基函数中的每一个对于具有零二阶和三阶导数。 X≥ξK
y=β0+β1x+β2([(x−k1)+3−(x−k4)+3k4−k1]−[(x−k3)+3−(x−k4)+3k4−k3])+β3([(x−k2)+3−(x−k4)+3k4−k2]−[(x−k3)+3−(x−k4)+3k4−k3])

我可以问一个更简单的问题吗：网站/书籍说：对于我的自然三次样条方法（即受限三次样条）w/4 节，我需要 4 个基函数。是 Beta_0、Beta_1*x 和 '4 more' 吗？或者确实只有 4 个 beta（从概念上讲，正如我上面所说的）？谢谢你。

感谢您的指导。我正在拟合样条曲线以及许多其他建模协变量，即解释变量。因此，仅用于样条的简单罐装包装是不够的。
当我拟合我的建模数据时，我期望的实际形状是一个非常积极的偏斜分布形式（我试图拟合具有拐点的奇特形状的示例数据，认为这将是对我的样条编码的一个很好的测试功能形式）。

我考虑过窃取函数形式和校准参数化（来自上面的 Python 或来自 R）——但它是三次样条，而不是自然三次样条。

我知道我的 ftn 如何与第一个结的 LHS 呈线性关系。我想下一步是让我看到各种术语取消，实际上我也会与最右边结的 RHS 成线性关系。还要寻找一个被称为“稳定”的 4 期基数。

3个回答

首先回答你的第二个问题：对于结，自然三次样条实际上只有自由度/基函数。KK

如果计算自由度，这很容易看出：个结将实线分成区间，这些区间中的分段三次函数具有个自由度。通过匹配每个结中的函数值、一阶和二阶导数，我们得到了个条件，因此个自由度。最后，“自然”样条条件表示我们希望样条在第一个和最后一个区间内是线性的，这会扼杀另一个自由度（这些区间中的二次和三次系数）。总而言之，我们剩下个自由度。KK+14(K+1)3KK+42+2K

对于您的第一个问题，如果不查看您的代码，真的很难判断出了什么问题。一般来说，我不建议您自己实现样条例程，除非您必须这样做；您使用的基础在数值上非常不稳定，并且通常代码有点挑剔才能正确。

但是，任何严肃的科学计算环境都内置了这些东西，因此通常您不需要自己编写。这是您在 Python/scipy 中的示例：

In [1]: %pylab
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]: n = 100

In [3]: x = linspace(0, 1, n)

In [10]: y = (((exp(1.2*x) + 1.5*sin(7*x))-1)/3) + normal(0, 0.15, size=n)

In [13]: from scipy.interpolate import LSQUnivariateSpline

In [15]: t = [0.2, 0.4, 0.6, 0.8]

In [16]: spl = LSQUnivariateSpline(x, y, t)

In [21]: plt.style.use('bmh')

In [22]: plot(x, spl(x), '-', x, y, 'o')

In [23]: plot(x, (((exp(1.2*x) + 1.5*sin(7*x))-1)/3), '--')

平滑 B 样条

蓝线是具有四个节点的平滑 B 样条，虚线是没有噪声的原始函数。

仅供参考，Hastie 等人的公式是正确的。我刚刚在 R 中自己实现了它，并将结果与​​我从 R 中的 splines::ns() 得到的结果进行了比较。基础并不相同，但它跨越了相同的空间，因为拟合值完全相同。

您的实现中某处存在错误。我采用了愚蠢的方法，为四个结实现了七个小函数：d1-d3、N1-N4，并且效果很好。

这是受限三次样条的实现（基于此matlab 代码）

import scipy.linalg as lin

def rcs(x,y,knots):
    n = len(y)
    k = knots
    X1 = x
    q = len(k)-1
    myX=np.zeros((n,len(knots)-2))

    for j in range(q-1):
        tmp1 = (x-k[j])**3 * (x>k[j])
    tmp2 = (x-k[q-1])**3 * (x>k[q-1])*(k[q]-k[j])
    XX= tmp1-tmp2/(k[q]-k[q-1])
        tmp1 = (x-k[q])**3 * (x>k[q])
        tmp2 = (k[q-1]-k[j])
    XX = XX+tmp1*tmp2/(k[q]-k[q-1])
    myX[:,j]=XX

    X = np.hstack( (np.ones((n,1)),np.reshape(X1,(n,1)),myX) )
    bhat = np.linalg.lstsq(X,y)[0]
    bhatt = np.zeros(len(knots)+1)
    bhatt[len(bhat)] = (bhat[2:]*(k[0:-2]-k[-1])).sum()
    bhatt[len(bhat)] = bhatt[len(bhat)] / (k[-1]-k[-2])
    bhatt = np.hstack([bhatt, 0])    
    bhatt[-1] = (bhat[2:]*(k[0:-2]-k[-2])).sum()
    bhatt[-1] = bhatt[-1] / (k[-2]-k[-1])
    bhat = np.hstack((bhat, bhatt[-2:]))
    return bhat

def speval(x,coefs,knots):
    tmp = coefs[0] + coefs[1]*x
    for k in range(len(knots)): tmp = tmp + coefs[k+2]*((x-knots[k])**3)*(x>knots[k])
    return tmp



import pandas as pd
x = np.random.randn(300)*np.sqrt(2)
e = np.random.randn(300)*np.sqrt(0.5)
y = np.sin(x)+e
df = pd.DataFrame([x,y]).T
df.columns = ['x','y']
df = df.sort_index(by='x')
print df.head()
knots=np.array([-5.5938, -3.7732, -1.9526, -0.1320, 1.6886, 3.5092, 5.3298]);
bhat = rcs(df.x,df.y,knots)
print bhat
df['spline'] = speval(df.x, bhat, knots)
df2 = df.set_index('x')
df2[['y','spline']].plot()
plt.hold(True)
for k in knots: plt.plot(k,speval(k,bhat,knots),'rd')

在此处输入图像描述

其它你可能感兴趣的问题
上一篇L-峰度的有效计算？ 下一篇二变量函数的自适应绘图z= f( x , y)z=f(x,y)算法伪代码？