本周，来自技术营销部的Richard Willey将会发表关于非参数拟合的客座博客。

内容

这周的博客是由一个常见的问题MATLAB中央．考虑以下问题:

你如何生成一个曲线拟合，最好地描述X和Y之间的关系?这个问题通常以两种不同的形式出现:

在本周的博客中，我将向你展示一种更好的方法来解决同样的问题——如何结合使用局部回归和交叉验证。

生成一些要使用的数据

为了使用更直观的方法，假设你有以下一组数据点:

s = RandStream(“mt19937ar”，“种子”, 1971);RandStream.setDefaultStream(年代);X = linspace(1,10,100);X = X';指定二阶傅里叶级数的参数W = .6067;A0 = 1.6345;A1 = -.6235;B1 = -1.3501;A2 = -1.1622;B2 = -.9443;%傅里叶2是X和Y之间的真实(未知)关系Y = a0 + a1 * cos (X * w) + b1 * sin (X * w) + a2 * cos (X 2 * * w) + b2 *罪(2 * X * w);添加一个噪声向量K = max(Y) - min(Y);噪声= Y + 2*K*randn(100,1);生成散点图散射(X,吵闹,“k”）;L2 =传说(“噪声数据样本”2);snapnow

非线性回归例子

如果你知道这个数据是用二阶傅里叶级数生成的，你可以使用非线性回归来建模Y = f(X)。

foo = fit(X，有噪声，“fourier2”）绘制结果持有在plot(foo) L3 =图例(“噪声数据样本”，非线性回归的2);持有从snapnow

foo =一般模型Fourier2: foo(x) = a0 + a1*cos(x*w) + b1*sin(x*w) + a2*cos(2*x*w) + b2*sin(2*x*w)系数(95%置信限):a0 = 1.734 (1.446, 2.021) a1 = -0.1998 (-1.065, 0.6655) b1 = -1.413 (-1.68， -1.146) a2 = -0.7688 (-1.752, 0.2142) b2 = -1.317 (-1.867， -0.7668) w = 0.6334 (0.5802, 0.6866)

非参数拟合实例

一切都很好，很简单……我们可以使用几行代码指定一个回归模型并绘制我们的曲线。我们的回归模型估计的系数与我们最初指定的系数非常接近。

问题来了:我们作弊了。当我们创建回归模型时，我们指定X和Y之间的真实关系应该使用二阶傅里叶级数来建模。如果我们无法获得这些信息呢?

下一段代码有点复杂，但是，它将只使用原始数据和一些CPU周期来解决相同的问题。

我们将使用一种称为LOWESS(局部散点图平滑)的平滑算法来生成我们的曲线拟合。LOWESS拟合的准确性取决于为平滑参数指定“正确的”值。我们将生成99个不同的LOWESS模型，使用0到1之间的平滑参数，看看哪个值生成最准确的模型。

这种蛮力技术经常遇到“过拟合”的问题。最小化“误差平方和”的LOWESS曲线适合数据集的信号和噪声成分。

为了防止过度拟合，我们将使用一种称为交叉验证的技术。我们要把数据集分成不同的训练集和测试集。我们将使用训练集中的数据生成预测模型，然后使用测试集中的数据测量模型的准确性。

Num = 99;span = linspace(.01，.99,num);Sse = 0(大小(跨度));Cp = cvpartition(100，“k”10);为j = 1:长度(跨度),f = @(火车、测试)规范(测试(:,2)——mylowess(火车、测试(:1),跨度(j))) ^ 2;sse(j) = sum(crossval(f，[X，噪声]，“分区”, cp));结束[minse,minj] = min(sse);Span = Span (minj);

mylowess是我写的一个函数，用来估计低平滑的值。我附上了一份mylowess.m供你使用。

类型mylowess.m

函数ys=mylowess(xy,xs,span) % mylowess低平滑，保留x值% ys=mylowess(xy,xs)返回x/y数据在%两列矩阵xy中的平滑版本，但在xs处计算平滑并返回在ys中的%平滑值。任何超出XY范围的值都被认为%等于最接近的值。如果nargin<3 ||为空(span) span = .3;对xy数据进行排序并获得平滑版本的xy = sortrows(xy);X1 = xy(:，1);Y1 = xy(:，2);Ys1 =光滑(x1,y1,span，'黄土');删除重复，这样我们就可以插值t = diff(x1)==0;x1 (t) = [];Ys1 (t) = []; % Interpolate to evaluate this at the xs values ys = interp1(x1,ys1,xs,'linear',NaN); % Some of the original points may have x values outside the range of the % resampled data. Those are now NaN because we could not interpolate them. % Replace NaN by the closest smoothed value. This amounts to extending the % smooth curve using a horizontal line. if any(isnan(ys)) ys(xsx1(end)) = ys1(end); end

评估我们的LOWESS拟合的准确性

接下来，让我们创建一个显示

原来的二阶傅里叶级数

这条曲线的噪声样本

非线性回归模型

LOWESS模型

情节(X, Y,“k”）;持有在散射(X,吵闹,“k”）;情节(foo,“r”) x = linspace(min(x)，max(x));线(x, mylowess (x,嘈杂的,x,跨度),“颜色”，“b”，“线型”，“- - -”，“线宽”2)传说(“干净数据”，噪声样本的，非线性回归的，“洛斯”2) hold从snapnow

正如您所看到的，LOWESS曲线很好地捕捉了系统的动态。我们能够生成一个LOWESS拟合，几乎和非线性回归一样准确，而不需要指定这个系统应该被建模为二阶傅里叶级数。

特别嘉宾出场由bootstrp

这是给那些坚持到最后的人的奖励。下一段代码使用“配对引导”生成LOWESS拟合的置信区间。

scatter(X, noise) f = @(xy) mylowess(xy,X,span);yboot2 = bootstrp(1000,f，[X，噪声])';mean黄土= mean(yboot2,2);h1 = line(X，均值黄土，“颜色”，“k”，“线型”，“- - -”，“线宽”2);std黄土= std(yboot2,0,2);h2 = line(X, mean黄土+2* std黄土，“颜色”，“r”，“线型”，“——”，“线宽”2);h3 = line(X, meanloess-2* std黄土，“颜色”，“r”，“线型”，“——”，“线宽”2);L5 =传说(“局部回归”，置信区间的2);snapnow

谁能提供一个很好的直观的解释如何引导工作?额外学分:

张贴你的答案在这里．

获取MATLAB代码

使用MATLAB®7.11发布