在 SciPy 中生成 B-Spline 基础，就像 R 中的 bs() [英] Generate a B-Spline basis in SciPy, like bs() in R

问题描述

对于 N 个一维数据 X，我想评估 K 个三次 B 样条上的每个点.在 R 中，有一个带有直观 API 的简单函数，称为

它不起作用，让我意识到我实际上并不知道这个函数在做什么.首先，它不会接受少于8个内部结，我不明白为什么.其次，它只认为样条是在 (1/3, 2/3)ish 范围内定义的，这可能意味着它出于某种原因忽略了前 3 个和最后 3 个结值?我需要打结吗?

任何帮助将不胜感激！

我已经解决了这个差异，实际上 BSpline 似乎忽略了结的前 3 个和最后 3 个值.我仍然很想知道为什么会出现这种差异，这样我就不会因为花在调试奇怪界面上的奇怪时间而感到难过.

为了后代，这里是产生基函数的代码

将 numpy 导入为 np导入 scipy.interpolate 作为 intrp导入 matplotlib.pyplot 作为 plt导入 patsy # 进行比较这些_knots = np.linspace(0,1,5)# 在 Patsy/R 中:善良而明智x = np.linspace(0., 1., 100)y = patsy.bs(x,knots=these_knots, degree=3)plt.subplot(1,2,1)plt.plot(x,y)plt.title('B样条基础')# 在 scipy: ?????numpyknots = np.concatenate(([0,0,0],these_knots,[1,1,1])) # 因为??y_py = np.zeros((x.shape[0], len(these_knots)+2))对于范围内的 i(len(these_knots)+2):y_py[:,i] = intrp.BSpline(numpyknots, (np.arange(len(these_knots)+2)==i).astype(float), 3, extrapolate=False)(x)plt.subplot(1,2,2)plt.plot(x,y_py)plt.title('在 SciPy')

解决方案

看起来您已经找到了答案，但要阐明为什么需要在边缘定义多个结，您可以阅读 scipy 文档.它们是使用 Cox-de Boor 递归公式定义的.该公式首先定义给定节点之间的相邻支持域，其常量值为 1(零阶).这些被卷积以获得更高阶的基函数.因此，两个域构成一个一阶基函数，三个域构成一个二阶基函数，四个域(= 5 个结点)构成一个在这 5 个结点范围内支持的三阶基函数.如果你想要 n 个 k = 3 的基函数，你需要有 (n+k+1) 个结点.

8 节的最小值使得 n >= k + 1，即 2 * (k+1).scipy 中的基区间 t[k] ... t[n] 是唯一可以定义全度基函数的范围.为了确保这个基本间隔到达外部结点，通常给两个末端结点的重数为 (k+1).scipy 可能只在您的其他"结果中显示了这个基本间隔.

请注意，您也可以使用

获得基函数

y_py[:,i] = intrp.BSpline.basis_element(numpyknots[i:i+5], extrapolate=False)(x)

这也消除了 x = 1 处的差异.

With N 1-dimensional data X, I would like to evaluate each point at K cubic B-splines. In R, there is a simple function with an intuitive API, called bs. There is actually a python package patsy which replicates this, but I can't use that package -- only scipy and such.

Having looked through the scipy.interpolate documentation on spline-related functions, the closest I can find is BSpline, or BSpline.basis_element, but how to get just the K basis functions is totally mysterious to me. I tried the following:

import numpy as np
import scipy.interpolate as intrp
import matplotlib.pyplot as plt
import patsy # for comparison

# in Patsy/R: nice and sensible
x = np.linspace(0., 1., 100)
y = patsy.bs(x, knots=np.linspace(0,1,4), degree=3)
plt.subplot(1,2,1)
plt.plot(x,y)
plt.title('B-spline basis')

# in scipy: ?????
y_py = np.zeros((x.shape[0], 6))
for i in range(6):
    y_py[:,i] = intrp.BSpline(np.linspace(0,1,10),(np.arange(6)==i).astype(float), 3, extrapolate=False)(x)


plt.subplot(1,2,2)
plt.plot(x,y_py)
plt.title('Something else')

It doesn't work, and makes me realise I don't actually know what this function is doing. First of all, it will not accept fewer than 8 interior knots, which I don't understand why. Secondly, it only thinks that the splines are defined within (1/3, 2/3)ish range, which maybe means that it is ignoring the first 3 and last 3 knot values for some reason? Do I need to pad the knots?

Any help would be appreciated!

EDIT: I have solved this discrepancy, indeed it seems like BSpline ignore the first 3 and last 3 values of knots. I'm still interested in knowing why there is this discrepancy, so that I feel less bad for the odd hour spent debugging a strange interface.

For posterity, here is the code that does produce the basis functions

import numpy as np
import scipy.interpolate as intrp
import matplotlib.pyplot as plt
import patsy # for comparison

these_knots = np.linspace(0,1,5)

# in Patsy/R: nice and sensible
x = np.linspace(0., 1., 100)
y = patsy.bs(x, knots=these_knots, degree=3)
plt.subplot(1,2,1)
plt.plot(x,y)
plt.title('B-spline basis')

# in scipy: ?????
numpyknots = np.concatenate(([0,0,0],these_knots,[1,1,1])) # because??
y_py = np.zeros((x.shape[0], len(these_knots)+2))
for i in range(len(these_knots)+2):
    y_py[:,i] = intrp.BSpline(numpyknots, (np.arange(len(these_knots)+2)==i).astype(float), 3, extrapolate=False)(x)

plt.subplot(1,2,2)
plt.plot(x,y_py)
plt.title('In SciPy')

解决方案

Looks like you already found the answer, but to clarify why these you need to define the multiple knots at the edges, you can read the scipy docs. They are defined using the Cox-de Boor recursive formula. This formula starts with defining neighbouring support domains between the given knot points with a constant value of 1 (zeroth order). These are convoluted to acquire the higher order basis functions. Hence two domains make one first order basis function, three domains make one second order basis function and four domains (= 5 knot points) make one third order basis function that is supported within the range of these 5 knot points. If you want n basis functions of degree k = 3, you will need to have (n+k+1) knot points.

The minimum of 8 knots is such that n >= k + 1, which gives 2 * (k+1). The base interval t[k] ... t[n] in scipy is the only range where you can define full degree basis functions. To make sure that this base interval reaches the outer knot points, the two end knots are usually given a multiplicity of (k+1). Probably scipy only showed this base interval in your 'Something else' result.

Note that you can also get the basis functions using

y_py[:,i] = intrp.BSpline.basis_element(numpyknots[i:i+5], extrapolate=False)(x)

this also removes the difference at x = 1.

这篇关于在 SciPy 中生成 B-Spline 基础，就像 R 中的 bs()的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！