Matlab无限接近0的数显示0,多变量灰色预测模型计算模型系数时，P矩阵奇异值接近零怎么解决？...

本帖最后由 王昕彤 于 2019-3-13 16:44 编辑

补充完整程序如下：

clc,clear all;

format long g;

fprintf(2,'OBGM(1,N)模型:\n');

X0=input('请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开),如:[2.874 3.278 3.307 3.390 3.679;7.04 7.645 8.075 8.53 8.774]\n:');

X00=input('请输入预测时的自变量序列(多个自变量序列之间用分号隔开),如:[911 12 13]\n:');

epsilon=0.5;

[n,m]=size(X0);

F=zeros(n,m);

for u=1:n

ago=0;

for v=1:m

ago=ago+X0(u,v);

F(u,v)=ago;

end

end

Y=X0';

Y=Y(:,1);

Y(1,:)=[];

B=F';

B(1,:)=[];B(:,1)=[];

Bm=[];Li=[];Co=[];

for k=1:m-1

Bm=[Bm,(epsilon-1)*F(1,k)-epsilon*F(1,k+1)];

Li=[Li,k+1];

Co=[Co,1];

end

B=[B,Bm'];

B=[B,Li'];

B=[B,Co'];

P=(B'*B)^(-1)*B'*Y;

len=length(P);

in_var_number=n-1;

a=P(in_var_number+1);

c=P(in_var_number+2);

d=P(in_var_number+3);

theta_1=1/(1+epsilon*a);

theta_2=1-a/(1+epsilon*a);

theta_3=c/(1+epsilon*a);

theta_4=d/(1+epsilon*a);

Simu1=[];

x01=X0(1,1);

Simu1=[Simu1,x01];

for k=2:m

s1=0;s2=0;s3=0;

for u=1:k-1

for i=2:n

s1=s1+theta_1*power(theta_2,u-1)*P(i-1)*F(i,k-u+1);

end

end

for v=0:k-2

s2=s2+power(theta_2,v)*((k-v)*theta_3+theta_4);

end

s3=power(theta_2,k-1)*x01;

Simu1=[Simu1,s1+s2+s3];

end

Xf=X0;

Xf(1,:)=[];

X0_F=[Xf,X00];

[g,h]=size(X0_F);

for u=1:g

ago=0;

for v=1:h

ago=ago+X0_F(u,v);

X0_F_1(u,v)=ago;

end

end

Fore1=[];

for k=m+1:h

s1=0;s2=0;s3=0;

for u=1:k-1

for i=2:n

s1=s1+theta_1*power(theta_2,u-1)*P(i-1)*X0_F_1(i-1,k-u+1);

end

end

for v=0:k-2

s2=s2+power(theta_2,v)*((k-v)*theta_3+theta_4);

end

s3=power(theta_2,k-1)*x01;

Fore1=[Fore1,s1+s2+s3];

end

Fore1=[F(1,m), Fore1];

Fore0=[];

for fo=1:length(Fore1)-1

ss= Fore1(fo+1)-Fore1(fo);

Fore0=[Fore0,roundn(ss,-3)];

end

A=zeros(n,5);

A(1,1)=1;

A(1,2)=X0(1,1);

A(1,3)=X0(1,1);

mp=0;

for k=2:length(X0)

A(k,1)=k;

A(k,2)=X0(1,k);

A(k,3)=roundn(Simu1(k)-Simu1(k-1),-3);

A(k,4)=roundn(A(k,3)-A(k,2),-3);

A(k,5)=roundn(100*abs(A(k,4))/A(k,2),-3);

mp=mp+100*abs(A(k,4))/A(k,2);

end

mp=roundn(mp/(length(X0)-1),-4);

disp('---------------------');

disp('(1)模型参数b2,b3,..,bn及a,c,d分别为:');

disp(roundn(P,-6));

disp('(2)误差检验表为:');

disp(' 序号 实际数据 模拟数据 残差 相对误差(%)');

disp(A);

disp('(3)平均相对模拟百分误差为(%):');

disp(mp);

str1='(4)未来[';

str2=int2str(length(Fore0));%数字转换为字符串

str3=']步的预测值分别为:';

str=[str1,str2,str3];%字符串连接

disp(str);

disp(Fore0);