多项式拟合正弦函数

%最小二乘法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变化情况');