下一个：，前一个：，上级：多项式运算 [目录][索引]

28.5 多项式插值

Octave 对各种插值提供了良好的支持，其中大部分在插值章节中进行了描述。上述章节中描述的函数的一个简单替代方案是将单个多项式或分段多项式（样条）拟合到给定的数据点。为了避免多项式剧烈波动，通常希望将低阶多项式拟合到数据中。这通常意味着需要在最小二乘意义下拟合多项式，而这正是polyfit函数所做的。

 
p = polyfit (x, y, n)
[p, s] = polyfit (x, y, n)
[p, s, mu] = polyfit (x, y, n)

返回使得点集 [x(:), y(:)] 的拟合最小二乘误差最小化的 n 次多项式 p(x) 的系数。

n 通常是一个 ≥ 0 的整数，用于指定近似多项式的阶数。如果 n 是一个逻辑向量，则将其用作掩码来选择性地强制使用或忽略相应的多项式系数。

多项式系数在行向量 p 中返回。输出 p 可以直接与 polyval 一起使用，以使用拟合的多项式来估计值。

可选的输出 s 是一个包含以下字段的结构体：

‘yf’

每个 x 值对应的多项式拟合值。

‘V’

用于计算多项式系数的范德蒙矩阵。

‘X（已弃用，将在 Octave 13 中移除）’

用于计算多项式系数的范德蒙矩阵。

‘R’

来自 QR 分解的三角因子 R。

‘C’

未缩放的协方差矩阵，形式上等于 v'*v 的逆，但以最小化舍入误差传播的方式计算。

‘df’

自由度。

‘normr’

残差的范数。

‘rsquared’

决定系数，即 R²。

第二个输出可供 polyval 用于计算预测值的统计误差限。特别地，p 系数的标准差由下式给出：

sqrt (diag (s.C)/s.df) * s.normr。

当存在第三个输出 mu 时，原始数据会被中心化和缩放，这可以提高拟合的数值稳定性。系数 p 与以下多项式关联：

xhat = (x - mu(1)) / mu(2)
其中 mu(1) = mean(x)，mu(2) = std(x)。

示例 1：逻辑 n 和整数 n 

f = @(x) x.^2 + 5;   # 数据生成函数
x = 0:5;
y = f (x);
## 将数据拟合到多项式 A*x^3 + B*x^1
p = polyfit (x, y, logical ([1, 0, 1, 0]))
⇒  p = [ 0.0680, 0, 4.2444, 0 ]
## 使用最高 x^3 的所有项将数据拟合到多项式
p = polyfit (x, y, 3)
⇒  p = [ -4.9608e-17, 1.0000e+00, -1.6906e-15, 5.0000e+00 ]

编程注意：如果期望的多项式阶数 n ≥ 数据点的数量，则解不唯一。此时，拟合计算 m = numel (x) - 1 个项。返回的多项式包含 x^n、x^n-1、…、x^n-m 以及常数项 x^0 的系数，其余系数均设为 0。

另请参阅： polyval、polyaffine、roots、vander、zscore。

当单个多项式不够好时，解决方案是使用多个多项式拼接在一起。splinefit 函数将分段多项式（样条）拟合到一组数据。

 
pp = splinefit (x, y, breaks)
pp = splinefit (x, y, p)
pp = splinefit (…, "periodic", periodic)
pp = splinefit (…, "robust", robust)
pp = splinefit (…, "beta", beta)
pp = splinefit (…, "order", order)
pp = splinefit (…, "constraints", constraints)

使用断点（节点）breaks 将分段三次样条拟合到带噪声的数据 x 和 y。

x 是一个向量，y 是一个向量或 N 维数组。如果 y 是 N 维数组，则 x(j) 与 y(:,…,:,j) 匹配。

p 是一个正整数，定义了沿 x 的区间数量，p+1 是断点的数量。每个区间中的点数相差不超过 1。

可选属性 periodic 是一个逻辑值，用于指定是否对样条应用周期性边界条件。周期的长度为 max (breaks) - min (breaks)。默认值为 false。

可选属性 robust 是一个逻辑值，用于指定是否应用稳健拟合以减少异常数据点的影响。将执行三次加权最小二乘迭代。权重根据之前的残差计算。异常值识别的灵敏度由属性 beta 控制。beta 的值限制在范围 0 < beta < 1 内。默认值为 beta = 1/2。接近 0 的值给予所有数据相等的权重。beta 的递增会减少异常数据的影响。接近 1 的值可能导致不稳定或秩亏。

拟合的样条以分段多项式 pp 的形式返回，可以使用 ppval 进行求值。

样条由阶数为 order 的多项式构成。默认是三次样条，order=3。具有 P 段的样条有 P+order 个自由度。在周期性边界条件下，自由度减少到 P。

可选属性 constraints 是一个结构体，用于指定拟合的线性约束。该结构体有三个字段："xc"、"yc" 和 "cc"。

"xc"

约束位置的 x 值向量。

"yc"

在位置 xc 处的约束值。默认是零数组。

"cc"

系数（矩阵）。默认是单位数组。行数限于分段多项式的阶数 order。

约束是阶数为 0 到 order-1 的导数的线性组合，形式如下： 

cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).

另请参阅： interp1、unmkpp、ppval、spline、pchip、ppder、ppint、ppjumps。

用于构造分段多项式的 breaks（或节点）数量是抑制输入数据 x 和 y 中噪声的重要因素。下面的示例演示了这一点。 

x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## 均匀断点
breaks = linspace (0, 2 * pi, 41); % 41 个断点，40 段
pp1 = splinefit (x, y, breaks);
## 从数据插值断点
pp2 = splinefit (x, y, 10);  % 11 个断点，10 段
## 绘图
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})

结果见图 图 28.1。

splinefit1

图 28.1：使用 41 个断点和 11 个断点拟合分段多项式的比较。断点数量多的拟合显示出底层函数中不存在快速波动。

由 splinefit 提供的分段多项式拟合具有直到 order-1 阶的连续导数。例如，三次拟合具有连续的一阶和二阶导数。以下代码演示了这一点： 

## 数据（200 个点）
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## 分段常数
pp1 = splinefit (x, y, 8, "order", 0);
## 分段线性
pp2 = splinefit (x, y, 8, "order", 1);
## 分段二次
pp3 = splinefit (x, y, 8, "order", 2);
## 分段三次
pp4 = splinefit (x, y, 8, "order", 3);
## 分段四次
pp5 = splinefit (x, y, 8, "order", 4);
## 绘图
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})

结果见图 图 28.2。

splinefit2

图 28.2：使用 8 个断点的分段常数、线性、二次、三次和四次多项式对带噪声数据拟合的比较。高阶解更准确地表示底层函数，但计算复杂度也更高。

当需要拟合的底层函数是周期性的时，splinefit 能够应用实现周期性拟合所需的边界条件。以下代码演示了这一点。 

## 数据（100 个点）
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## 无约束
pp1 = splinefit (x, y, 10, "order", 5);
## 周期性边界
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## 绘图
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})

结果见图 图 28.3。

splinefit3

图 28.3：有无周期性边界条件下分段多项式拟合带噪声周期函数的比较。

还可以添加更复杂的约束。例如，下面的代码展示了一个周期拟合，其端点的值被钳制，以及第二个端点处铰接的周期拟合。 

## 数据（200 个点）
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## 断点
breaks = linspace (0, 2 * pi, 10);
## 端点钳制，y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## 铰接式周期端点，y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## 绘图
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})

结果见图 图 28.4。

splinefit4

图 28.4：两个周期分段三次拟合对带噪声周期信号的比较。一个拟合的端点被钳制，第二个拟合的端点被铰接。

splinefit 函数还提供了便捷的稳健（robust）拟合功能，可以减少异常数据的影响。在下面的示例中，提供了三种不同的拟合，两种具有不同程度的异常值抑制，第三种展示了非稳健解。 

## 数据
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## 添加异常值
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## 使用铰接条件拟合样条
con = struct ("xc", [0, 2*pi]);
## 稳健拟合，beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## 稳健拟合，beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## 无稳健拟合
pp3 = splinefit (x, y, 8, "constraints", con);
## 绘图
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
         "robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto

结果见图 图 28.5。

splinefit6

图 28.5：两种不同稳健拟合水平（beta = 0.25 和 0.75）对带噪声数据与异常值数据的比较。还包含一个没有稳健拟合的常规拟合（beta = 0）。

一种非常特殊的多项式插值形式是 Padé 逼近。对于控制系统，连续时间延迟可以非常简单地用逼近器建模。

 
[num, den] = padecoef (T)
[num, den] = padecoef (T, N)

计算连续时间延迟 T 的 N 阶 Padé 逼近，以传递函数形式表示。exp (-sT) 的 Padé 逼近由以下方程定义： 

             Pn(s)
exp (-sT) ~ -------
              Qn(s)

其中 Pn(s) 和 Qn(s) 都是由以下表达式定义的 N 阶有理函数： 

         N    (2N - k)!N!        k
Pn(s) = SUM --------------- (-sT)
        k=0 (2N)!k!(N - k)!

Qn(s) = Pn(-s)

输入 T 和 N 必须是非负数值标量。如果未指定 N，则默认为 1。

输出行向量 num 和 den 分别包含分子和分母系数，按 s 的降幂排列。两者都是 N 阶多项式。

例如： 

t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒  num =

      1.0000e-04  -2.0000e-02   1.8000e+00  -8.4000e+01   1.6800e+03

⇒  den =

      1.0000e-04   2.0000e-02   1.8000e+00   8.4000e+01   1.6800e+03

ppval 函数用于求值由 mkpp 或其他方式创建的分段多项式，而 unmkpp 返回关于分段多项式的详细信息。

以下示例展示了如何将两个线性函数和一个二次函数组合成一个函数。每个函数在相邻的区间上表示。 

x = [-2, -1, 1, 2];
p = [ 0,  1, 0;
      1, -2, 1;
      0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);
 
pp = mkpp (breaks, coefs)
pp = mkpp (breaks, coefs, d)

根据样本点 breaks 和系数 coefs 构造分段多项式 (pp) 结构体。

breaks 必须是严格递增值的向量。区间数由 ni = length (breaks) - 1 给出。

当 m 是多项式阶数时，coefs 的大小必须为：ni-by-(m + 1)。

coefs 的第 i 行 coefs(i,:) 包含第 i 个区间上多项式的系数，从最高阶 (m) 到最低阶 (0) 排列。

coefs 也可以是一个多维数组，用于指定向量值或数组值多项式。在这种情况下，多项式阶数 m 由 coefs 最后一个维度的长度定义。第一维的大小由标量或向量 d 给出。如果未指定 d，则默认为 1。在这种情况下，p(r, i, :) 包含定义在区间 i 上的第 r 个多项式的系数。在任何情况下，coefs 都会被重塑为大小为 [ni*prod(d) m] 的二维矩阵。

编程注意：ppval 在 xi - breaks(i) 处求值多项式，即它从 xi 中减去当前区间的下端点。在使用 mkpp 创建分段多项式对象时，必须考虑到这一点。

另请参阅： unmkpp、ppval、spline、pchip、ppder、ppint、ppjumps。

 
[x, p, n, k, d] = unmkpp (pp)

提取分段多项式结构体 pp 的组成部分。

此函数是 mkpp 的逆操作：它提取创建分段多项式结构体 pp 所需的 mkpp 输入。以下代码使这种关系更加明确： 

[breaks, coefs, numinter, order, dim] = unmkpp (pp);
pp2  = mkpp (breaks, coefs, dim);

以这种方式获得的分段多项式结构体 pp2 与原始 pp 相同。也可以通过直接访问结构体 pp 的字段来获得相同的结果。

各组成部分如下：

x

样本点或断点。

p

样本区间中点的多项式系数。p(i, :) 包含区间 i 上从高次到低次排列的多项式系数。如果 d > 1，则 p 是一个大小为 [n*prod(d) m] 的矩阵，其中 i + (1:d) 行是区间 i 中所有 d 个多项式的系数。

n

多项式段或区间的数量，n = length (x) - 1。

k

多项式的阶数加 1。

d

每个区间上定义的多项式数量。

另请参阅： mkpp、ppval、spline、pchip。

 
yi = ppval (pp, xi)

在点 xi 处求值分段多项式结构体 pp。

如果 pp 描述的是一个标量多项式函数，则结果是与 xi 形状相同的数组。否则，如果 xi 是向量，结果的大小为 [pp.dim, length(xi)]；如果 xi 是多维数组，则结果为 [pp.dim, size(xi)]。

另请参阅： mkpp、unmkpp、spline、pchip。

 
ppd = ppder (pp)
ppd = ppder (pp, m)

计算分段多项式结构体 pp 的 m 阶分段导数。

如果省略 m，则计算一阶导数。

另请参阅： mkpp、ppval、ppint。

 
ppi = ppint (pp)
ppi = ppint (pp, c)

计算分段多项式结构体 pp 的积分。

c（如果给定）是积分常数。

另请参阅： mkpp、ppval、ppder。

 
jumps = ppjumps (pp)

求值分段多项式的边界跳跃。

如果有 n 个区间，且 pp 的维度为 d，则结果数组的维度为 [d, n-1]。

另请参阅： mkpp。

下一个：其他函数，前一个：导数/积分/变换，上级：多项式运算 [目录][索引]

版权所有 © 2024-2026 Octave中文网

ICP备案/许可证号：黑ICP备2024030411号-2