图解机器学习第二部分第六章——鲁棒学习Huber损失最小化

n =10;%50 个节点
N=1000;%1000 个节点
x=linspace(-3,3,n)';%在 -3 到 3 上找 50 个节点
X=linspace(-4,4,N)';%在 -3 到 3 上找 1000 个节点 去拟合这 50 个节点
y=x+0.2*randn(n,1);%通过函数f(x) 得到 y。
y(n)=-4;% 设置一个异常值

p(:,1)=ones(n,1); %线性模型 a + bx
p(:,2)=x;
t0=p\y; %给参数设置一个适当的初值
e=1; %将阈值设置为 1
for o=1:1000
    r=abs(p*t0-y);% 计算残差的绝对值
    w=ones(n,1); %计算权重矩阵
    w(r>e)=e./r(r>e);
    t=(p'*(repmat(w,1,2).*p))\(p'*(w.*y));
    if norm(t-t0)<0.001,break,end
    t0=t;
end
P(:,1)=ones(N,1);
P(:,2)=X;
F=P*t;

figure(1);
clf;
hold on;
axis([-4 4 -4.5 3.5]);% X 轴的范围是-2.8 ~ 2.8。Y 轴的范围是 -0.5 ~ 1.2
plot(X,F,'g-');
plot(x,y,'bo');

这段代码主要是实现了一个基于最小二乘法的线性回归拟合过程，并使用散点图和拟合曲线进行可视化。

首先，代码定义了两个变量 n 和 N 分别表示节点的数量，其中 n = 10 表示有 10 个节点，N = 1000 表示有 1000 个节点。

接下来，代码使用 linspace 函数生成了两个向量 x 和 X，分别代表节点的位置。其中 x 通过在 -3 到 3 的范围内均匀生成 10 个节点，而 X 则生成了 1000 个节点，范围是 -4 到 4。

然后，代码通过 y = x + 0.2 * randn(n,1) 计算了对应节点位置 x 的函数值 y，并添加了一个正态分布噪声（均值为 0，标准差为 0.2）。最后，代码将第 n 个节点的函数值设置为 -4，模拟了一个异常值。

接下来的代码段开始进行回归拟合过程：

p(:,1) = ones(n,1) 和 p(:,2) = x 定义了一个矩阵 p，其中第一列是全为 1 的向量，第二列是节点位置 x。这里的 p 用于表示线性模型 a + bx 中的系数。
t0 = p\y 使用最小二乘法求解线性模型的参数，即拟合直线的截距和斜率的初值。
e = 1 设置了一个阈值，用于确定权重矩阵中的截断点。
接下来的循环执行拟合过程。在每次迭代中：
r = abs(p * t0 - y) 计算了残差的绝对值，即模型拟合值与实际值之间的差异。
w = ones(n,1) 初始化权重矩阵为全 1。
w(r > e) = e ./ r(r > e) 对于残差大于阈值 e 的节点，将对应的权重设为 e 除以残差。
t = (p' * (repmat(w,1,2) .* p)) \ (p' * (w .* y)) 使用加权最小二乘法更新模型参数。
if norm(t - t0) < 0.001, break, end 判断当前参数与上一次参数之间的差异是否小于阈值，如果是，则终止循环。
t0 = t 更新参数的值。

最后，代码定义了矩阵 P，用于表示拟合曲线的节点位置。与之前的 p 不同，P 是基于 X 生成的。
F = P * t 计算了拟合曲线的函数值。
最后的可视化部分使用 plot 函数绘制了拟合曲线（绿色实线）和原始数据点（蓝色圆点）。