優化GM(1,1)

主要的改進手段：優化微分方程的初始條件

data = csvread('D:\matlab2016a\demo\data\air.csv');
num = size(data,1);
predictnum = 6;
train = data(1:num-predictnum,:);
test = data((num-predictnum+1):end,:);

Y = train(2:end,:);

for i = 1:size(train,1)
    X(i,:) = sum(train(1:i,:));
end

for j = 2:size(X,1)
    Z(j-1,:) = (X(j,:)+X(j-1,:))/2;
end
B = [-Z,ones(size(Z,1),1)];

ab = inv(B'*B)*B'*Y;
a = ab(1,:);
b = ab(2,:);

for i = 2:num
    predict(i-1,:) = (X(end,:)-b/a)*(1-exp(a))*exp(-a*(i-size(train,1)));
end

predict_train = predict(1:size(train,1)-1,:);
predict_test = predict((size(train,1)):end,:);

train_rmse = sqrt(sum((predict_train-train(2:end,:)).^2)/length(predict_train));
train_mape = sum(abs((predict_train-train(2:end,:))./train(2:end,:)))/length(predict_train);
test_rmse = sqrt(sum((predict_test-test).^2)/length(predict_test));
test_mape = sum(abs((predict_test-test)./test))/length(predict_test);

disp(['Train rmse = ',num2str(train_rmse),' Train mape = ',num2str(train_mape*100)])
disp(['Test rmse = ',num2str(test_rmse),' Test mape = ',num2str(test_mape*100)])

(二)OICGM(1,1)

data = csvread('D:\matlab2016a\demo\data\air.csv');
num = size(data,1);
predictnum = 6;
train = data(1:num-predictnum,:);
test = data((num-predictnum+1):end,:);

Y = train(2:end,:);

for i = 1:size(train,1)
    X(i,:) = sum(train(1:i,:));
end

for j = 2:size(X,1)
    Z(j-1,:) = (X(j,:)+X(j-1,:))/2;
end
B = [-Z,ones(size(Z,1),1)];

ab = inv(B'*B)*B'*Y;
a = ab(1,:);
b = ab(2,:);

count = 1:size(train,1);
weights = count./sum(count');
s = (weights*X-b/a)*(1-exp(a));

r = 0;
for i = 1:size(train,1)
    r = r + exp(-a*i)*train(i,:);
end

t = 0;
for i = 1:size(train,1)
    t = t + exp(-2*a*i);
end

beta = log(r/(s*t))/a;

for i = 2:num
    predict(i-1,:) = (weights*X-b/a)*(1-exp(a))*exp(-a*(i-beta));
end

predict_train = predict(1:size(train,1)-1,:);
predict_test = predict((size(train,1)):end,:);

train_rmse = sqrt(sum((predict_train-train(2:end,:)).^2)/length(predict_train));
train_mape = sum(abs((predict_train-train(2:end,:))./train(2:end,:)))/length(predict_train);
test_rmse = sqrt(sum((predict_test-test).^2)/length(predict_test));
test_mape = sum(abs((predict_test-test)./test))/length(predict_test);

disp(['Train rmse = ',num2str(train_rmse),' Train mape = ',num2str(train_mape*100)])
disp(['Test rmse = ',num2str(test_rmse),' Test mape = ',num2str(test_mape*100)])