matlab算法记录学习

一、 残差自相关性分析，以及如何消除残差

详情建模数第六章

clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3;160.7;164.2;165.6;168.7;171.7]; 
rstool(x,y)


%数据Dw检验，相关性检验，函数考虑时间，证实时间是否有滞后性
a=ones(numel(x),1);%有多少个数据，要改
b=0.176*x-1.455*a;%预测函数值，要改
c=y-b   % 残差=实际值-预测值
c1=c(2:end,1);%残差e(t)
c2=c(1:end-1,1);%e(t-1)
plot(c1,c2,'*')%残差散点图
xlabel('e(i-1)'),
ylabel('e(i)')
hold on
%画横纵坐标
d=0;
d1=-0.15:0.001:0.25;
plot(d,d1,'.r')
hold on
plot(d1,d,'.r')

% DW的程序：
num=c(1:19)'*c(2:20);
den=sum(c(2:20).^2);
p = num/den%p值
DW=2*(1-p)%DW<dl,残差存在明显自相关,接着消除随机偏差项自相关后的回归模型。

clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3;160.7;164.2;165.6;168.7;171.7]; 
rstool(x,y)


%数据Dw检验，相关性检验，函数考虑时间，证实时间是否有滞后性
a=ones(numel(x),1);%有多少个数据，要改
b=0.176*x-1.455*a;%预测函数值，要改
c=y-b   % 残差=实际值-预测值
c1=c(2:end,1);%残差e(t)
c2=c(1:end-1,1);%e(t-1)
plot(c1,c2,'*')%残差散点图
xlabel('e(i-1)'),
ylabel('e(i)')
hold on
%画横纵坐标
d=0;
d1=-0.15:0.001:0.25;
plot(d,d1,'.r')
hold on
plot(d1,d,'.r')

% DW的程序：
num=c(1:19)'*c(2:20);
den=sum(c(2:20).^2);
p = num/den%p值
DW=2*(1-p)%DW<dl,残差存在明显自相关,接着消除随机偏差项自相关后的回归模型。

% DW模型优化后的程序：
clc,clear,close all
y=[20.96;21.40;21.96;21.52;22.39;22.76;23.48;23.66;24.10;24.01;24.54;24.30;25.00;25.64;26.36;26.98;27.52;27.78;28.24;28.78];
y1=y(2:length(y),1);
y2=y(1:length(y)-1,1);
y3=y1-0.71133*y2;%改y值,,0.71133为p值

x=[127.3;130.0;132.7;129.4;135.0;137.1;141.2;142.8;145.5;145.3;148.3;146.4;150.2;153.1;157.3; 160.7;164.2;165.6;168.7;171.7];
x1=x(2:length(x),1);
x2=x(1:length(x)-1,1);
x3=x1-0.71133*x2;

a=ones(length(x)-1,1);
y4=y3-(0.1763*x3-0.42*a)  % 残差,拟合得式子
y5=y4(2:length(y4),1);%残差e(t)
y6=y4(1:length(y4)-1,1);%残差e(t-1)
y7=y5-y6; %差值

y8=sum(y7.^2);
y9=sum(y5.^2);
DW=y8/y9
p=1-DW/2

二、时间序列线性二次平均预测法

%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%% 数据处理
y=y';
for i=1:numel(y)-2
    m1(i+2)=(y(i)+y(i+1)+y(i+2))/3;%一次移动平均
end
for i=3:numel(y)-2
    m2(i+2)=(m1(i)+m1(i+1)+m1(i+2))/3;%二次移动平均
end
m1=m1';
m2=m2';
a=2*m1-m2;
b=m1-m2;

%预测方程y(t+T)=at+bt*T,看书532页
at=a(length(y),1)%
bt=b(length(y),1)
for T=1:6
    y1(T,:)=at+bt*T;%预测将来数后六位
end
%% 预测值
y1
m2=[m2;y1];
for n=5:length(y)
  et(n,:)=abs(((m2(n)-y(n))/y(n))*100);%相对偏差分析
end
%% 画图
t1=5:1:length(x);
t2=5:1:length(x)+6;%预测后六年
y=y(5:t1(end));
m2=m2(5:t2(end));
plot(t1,y,'o',t2,m2)  %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

clc,clear,close all;
x1=1:1:39;%输入源数据
for i=1:numel(x1)-2
    m1(i+2)=(x1(i)+x1(i+1)+x1(i+2))/3;%一次移动平均
end
for i=3:numel(x1)-2
    m2(i+2)=(m1(i)+m1(i+1)+m1(i+2))/3;%二次移动平均
end
m1=m1';
m2=m2';
a=2*m1-m2;
b=m1-m2;

%预测方程y(t+T)=at+bt*T,看书532页
at=a(39,1)%
bt=b(39,1)
for T=1:6
    y(:,T)=at+bt*T;%预测将来数后六位
end
y

三、 灰度预测模型

%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%% 处理
syms a b;
c=[a b]';
A=y';%源数据
B=cumsum(A);  % 原始数据累加
n=length(A);
for i=1:(n-1)
    C(i)=(B(i)+B(i+1))/2;  % 生成累加矩阵
end
% 计算待定参数的值
D=A;D(1)=[];
D=D';
E=[-C;ones(1,n-1)];
c=inv(E*E')*E*D;
c=c';
a=c(1);b=c(2);

% 预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
    F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a ;
end
G=[];G(1)=A(1);

for i=2:(n+10)
    G(i)=F(i)-F(i-1); %获得预测出来的数据，预测十年
end 
t1=x';
t2=t1(1):t1(end)+10;
G, a, b % 输出预测值，发展系数和灰色做用量
%% 画图
plot(t1,A,'o',t2,G)  %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

4 一次指数平滑

%读取数据
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%一次平滑预测法
n=length(y);
alpha=[0.1,0.3,0.9];%长度
m=length(alpha);
yhat(1,1:m)=(y(1)+y(2))/2;
for i=2:n
    yhat(i,:)=alpha*y(i-1)+(1-alpha).*yhat(i-1,:);%一次指数平滑预测法
end
yhat;%预测值
et=(repmat(y,1,m)-yhat);%偏差
err=sqrt(mean((repmat(y,1,m)-yhat).^2))%预测偏差总和
yhatnaxt=alpha*y(n)+(1-alpha).*yhat(n,:)%预测下一时刻值
%画图
t=x';
plot(t,y,'o',t,yhat(:,3))  %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

1. 二次指数平滑

%% 二次平滑预测法
clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%数据处理
yx =y;
yt=yx;
n=length(yt);
alpha=0.8;
st1(1)=yt(1); 
st2(1)=yt(1);
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
end
a=2*st1-st2;
b=alpha/(1-alpha)*(st1-st2);
yhat=a+b;
yhat=yhat';
err=sqrt(mean((y-yhat).^2));%预测偏差总和
%画图
t=x';
plot(t,y,'o',t,yhat)  %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

1. 三次指数平滑

%% 三次指数平滑
clc,clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
yx =y;
yt=yx;
n=length(yt);
alpha=0.7;
st1_0=mean(yt(1:3));
st2_0=st1_0; 
st3_0=st1_0;
st1(1)=alpha*yt(1)+(1-alpha)*st1_0;
st2(1)=alpha*st1(1)+(1-alpha)*st2_0;
st3(1)=alpha*st2(1)+(1-alpha)*st3_0;
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
    st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1);
end
st1=[st1_0,st1]; 
st2=[st2_0,st2];
st3=[st3_0,st3];
a=3*st1-3*st2+st3;
b=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3);
c=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3);
yhat=a+b+c;
yhat=yhat';
err=sqrt(mean((y-yhat(1:end-1)).^2));%预测偏差总和
%画图
t=x';
plot(t,y,'o',t,yhat(1:end-1))  %原始数据与预测数据的比较
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

5画图线性拟合

clc;clear;
inputfile = '../MATLAB/1.xlsx' ; % 销量数据文件
outputfile ='../MATLAB/2.xlsx';  % 插值后数据存放
% 读入数据
[num,txt,raw] = xlsread(inputfile);
xname=txt(:,1);%时间名
x = num(:,1);
yname=txt(:,2);%数值名
y=num(:,2);
ti=txt(1,3);%标题
%%
cftool
%% 画图
plot(fittedmodel,x,y,'o');
xlabel(xname);
ylabel(yname);
title(ti);
set (gcf,'Position',[400,100,500,300], 'color','w')

6灰色关联度模型

CPI与什么因数有关

%%灰色关联度模型
clc,clear,close all
inputfile = '../MATLAB/hsgl.xlsx' ; % 销量数据文件
[num,txt,raw] = xlsread(inputfile);

r=size(num);
y=num(2:r,:);%源数据
y1=mean(y');
y1=y1';
[row,col]=size(y);%获取矩阵的行数，列数
for i=1:row
    for j=1:col
        y2(i,j)=y(i,j)/y1(i);%初值像矩阵
    end
end
for i=2:row
    for j=1:col
        y3(i-1,j)=abs(y2(i,j)-y2((i-1) ,j));%差序列
    end
end
a=1;b=0;
for i=1:row-1
       for j=1:col
            if(y3(i,j)<=a)  a=y3(i,j);%计算最大值
                elseif(y3(i,j)>=b) b=y3(i,j);%计算最小值
            end
        end
end
for i=1:row-1
    for j=1:col
            y4(i,j)=(a+0.5*b)/(y3(i,j)+0.5*b);%计算关联系数
    end
end
y5=sum(y4')/(col-1);%灰色关联度
y5=y5'