风速威布尔分布和ARMA预测模型matlab程序

clc
clear
%% 1.计算风速weibull分布
% 数据处理
load data;
mu=mean(speed);%原始数据的统计参数
sigma=sqrt(var(speed));
% 计算威布尔分布参数
parmhat=wblfit(speed);
k=parmhat(2);
c=parmhat(1);
% k=(sigma/mu)^-1.086;
% c=mu/gamma(1+1/k);
% 威布尔分布拟合
[y,x]=hist(speed,ceil(max(speed)/0.5));%x是区间中心数,组距-1.5
prob1=y/8760/0.5;%计算原始数据概率密度，频数除以数据种数，除以组距prob2=(k/c)*(x/c).^(k-1).*exp(-(x/c).^k);%威布尔分布
figure(1)
title('Weibull分布拟合图');
bar(x,prob1,1)
hold on
plot(x,prob2,'r')
legend('历史数据','Weibull拟合结果')
% legend('Weibull拟合结果')
hold off
save('result_weibull.mat')
%% 2.ARMA模型预测风速
clc
clear
load data
y=speed(1:300);
Data=y; %共300个数据
SourceData=Data(1:250,1); %前250个训练集
step=50; %后50个测试
TempData=SourceData;
TempData=detrend(TempData);%去趋势线
TrendData=SourceData-TempData;%趋势函数
%--------差分，平稳化时间序列---------
H=adftest(TempData);
difftime=0;
SaveDiffData=[];
while ~H
SaveDiffData=[SaveDiffData,TempData(1,1)];
TempData=diff(TempData);%差分，平稳化时间序列
difftime=difftime+1;%差分次数
H=adftest(TempData);%adf检验，判断时间序列是否平稳化
end
%---------模型定阶或识别--------------
u = iddata(TempData);
test = [];
for p = 1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q = 1:5 %移动平均对应ACF
m = armax(u,[p q]);
AIC = aic(m); %armax(p,q),计算AIC
test = [test;p q AIC];
end
end
for k = 1:size(test,1)
if test(k,3) == min(test(:,3)) %选择AIC值最小的模型
p_test = test(k,1);
q_test = test(k,2);
break;
end
end
%------1阶预测-----------------
TempData=[TempData;zeros(step,1)];
n=iddata(TempData);
%m = armax(u(1:ls),[p_test q_test]); %armax(p,q),[p_test q_test]对应AIC值最小，自动回归滑动平均模型
m = armax(u,[p_test q_test]);
% -------------------------------------------
P1=predict(m,n,1);
PreR=P1.OutputData;
PreR=PreR';
Noise.std=sqrt(m.NoiseVariance);
e=normrnd(0,Noise.std,1,300);
for i=251:300
PreR(i)=-m.A(2:p_test+1)*PreR(i-1:-1:i-p_test)'+m.C(1:q_test+1)*e(i:-1:i-q_test)'; end
% -------------------------------------------
%----------还原差分-----------------
if size(SaveDiffData,2)~=0
for index=size(SaveDiffData,2):-1:1
PreR=cumsum([SaveDiffData(index),PreR]);
end
end
%-------------------预测趋势并返回结果----------------
mp1=polyfit([1:size(TrendData',2)],TrendData',1);
xt=[];
for j=1:step
xt=[xt,size(TrendData',2)+j];
end
TrendResult=polyval(mp1,xt);
PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2));
tempx=[TrendData',TrendResult]+PreR; % tempx为预测结果
plot(tempx,'r-.');
hold on
plot(Data,'b');
legend('ARMA拟合时序曲线','实际时序风速');
save('resultarma.mat');。