%最小二乘法2.0
%分为带正则项和不带正则项。
close all;
clear all;
N=9;%最大特征值数
T = pi/2; %周期
step = T/100; %采样步长=T/10=pi/20;
lamda=0.1;%λ参数
t = (0 : step : 2*T);%取样t的函数值 采样频率
y = sin(T*t); %取样y的值，产生周期为pi/2的sin函数
%%
%显示图片，分别是原图，叠加高斯白噪声的正弦波取的点集。
figure(1);
hold on
plot(t,y,'b');
z1=0.35*randn(1,length(t)); %产生方差 N（0，0.12）高斯白噪声 （b=0.01/0.1/1
y2=y+z1; %叠加高斯白噪声的正弦波
plot(t,y2,'ro');%展示加过噪声后的图像。
legend('原图','叠加高斯白噪声的正弦波');
%其中t时自变量值，y2是生成的采集点。
%取测试样本
%%
%这里的机器学习结果存在x0,y0里面
%加过噪声的存在t,y2里
%原来的结果存在t,y里也就是测试样本
figure(3);
[~,k]=size(t);
w=zeros(1,N);
for n=1:N
X0=zeros(n+1,k);
for k0=1:k %构造矩阵X0
for n0=1:n+1
X0(n0,k0)=t(k0)^(n+1-n0);
end
end
X=X0';
ANSS=(X'*X)\X'*y2';%用函数解方程求逆了……
for i=1:n+1 %answer矩阵存储每次求得的方程系数，按列存储
answer(i,n)=ANSS(i);
end
x0=0 : step : 2*T;
y0=ANSS(1)*x0.^n ;%根据求得的系数初始化并构造多项式方程，求解loss
for num=2:1:n+1
y0=y0+ANSS(num)*x0.^(n+1-num);
end
%根据求得的系数初始化并构造多项式方程，在下面图里显示出来，
%为了保证平滑，取点更多
x0draw=0 : step*0.1 : 2*T;
y0draw=ANSS(1)*x0draw.^n;
for num=2:1:n+1
y0draw=y0draw+ANSS(num)*x0draw.^(n+1-num);
end
subplot(3,3,n);
hold on
plot(x0draw,y0draw,'r');
plot(t,y,'g');
%下面计算误差,存在w里
w(1,n)=1/(2*k)*sum((y-y0).^2);
end
figure(2);
plot(x0draw,y0draw,'r')
hold on
plot(t,y,'g');
figure(4)
plot(1:1:N,w,'r');
%%
%加入正则项
figure(3);
[~,k]=size(t);
wz=zeros(1,N);
for n=1:N
X0z=zeros(n+1,k);
for k0=1:k %构造矩阵X0
for n0=1:n+1
X0z(n0,k0)=t(k0)^(n+1-n0);
end
end
Xz=X0z';
ANSSZ=(Xz'*Xz+lamda.*eye(n+1))\Xz'*y2';%解方程求逆了……
for i=1:n+1 %answer矩阵存储每次求得的方程系数，按列存储
answerz(i,n)=ANSSZ(i);
end
x0z=0 : step : 2*T;
y0z=ANSSZ(1)*x0z.^n ;%根据求得的系数初始化并构造多项式方程，求解loss
for num=2:1:n+1
y0z=y0z+ANSSZ(num)*x0z.^(n+1-num);
end
%根据求得的系数初始化并构造多项式方程，在下面图里显示出来，
%为了保证平滑，取点更多
x0drawz=0 : step*0.1 : 2*T;
y0drawz=ANSSZ(1)*x0drawz.^n;
for num=2:1:n+1
y0drawz=y0drawz+ANSSZ(num)*x0drawz.^(n+1-num);
end
subplot(3,3,n)
hold on
plot(x0drawz,y0drawz,'b')
plot(t,y2,'ro');%展示加过噪声后的图像。
legend('不带正则项拟合曲线','原曲线','带正则项拟合曲线','所取加过噪声的点');
%下面计算误差,存在wz里
wz(1,n)=1/(2*k)*sum((y-y0z).^2);
end
figure(2);
plot(x0drawz,y0drawz,'b');
plot(t,y2,'ro');%展示加过噪声后的图像。
legend('没有加正则的9阶','原图','加入正则的9阶','所取的点');
wz
figure(4)
hold on
plot(1:1:N,wz,'b*')
legend('未加正则的loss变化情况','加正则的loss变化情况');
