matlab算法記錄學習

時間 2019-12-05

原文原文鏈接

一、殘差自相關性分析，以及如何消除殘差

詳情建模數第六章函數

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')

二次指數平滑優化

%% 二次平滑預測法
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')

三次指數平滑spa

%% 三次指數平滑
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與什麼因數有關code

%%灰色關聯度模型
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'