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

clcclear%%1.计算风速weibull分布%数据处理loaddata;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.5prob1=y/8760/0.5;%计算原始数据概率密度，频数除以数据种数，除以组距prob2=(k/c)*(x/c).^(k-1).*exp(-(x/c).^k);%威布尔分布figure(1)title('Weibull分布拟合图');bar(x,prob1,1)holdonplot(x,prob2,'r')legend('历史数据','Weibull拟合结果')%legend('Weibull拟合结果')holdoffsave('result_weibull.mat')%%2.ARMA模型预测风速clcclearloaddatay=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~HSaveDiffData=[SaveDiffData,TempData(1,1)];TempData=diff(TempData);%差分，平稳化时间序列difftime=difftime+1;%差分次数H=adftest(TempData);%adf检验，判断时间序列是否平稳化end%---------模型定阶或识别--------------u=iddata(TempData);test=[];forp=1:5%自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12forq=1:5%移动平均对应ACFm=armax(u,[pq]);AIC=aic(m);%armax(p,q),计算AICtest=[test;pqAIC];endendfork=1:size(test,1)iftest(k,3)==min(test(:,3))%选择AIC值最小的模型p_test=test(k,1);q_test=test(k,2);break;endend%------1阶预测-----------------TempData=[TempData;zeros(step,1)];n=iddata(TempData);%m=armax(u(1:ls),[p_testq_test]);%armax(p,q),[p_testq_test]对应AIC值最小，自动回归滑动平均模型m=armax(u,[p_testq_test]);%-------------------------------------------P1=predict(m,n,1);PreR=P1.OutputData;PreR=PreR';Noise.std=sqrt(m.NoiseVariance);e=normrnd(0,Noise.std,1,300);fori=251:300PreR(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%-------------------------------------------%----------还原差分-----------------ifsize(SaveDiffData,2)~=0forindex=size(SaveDiffData,2):-1:1PreR=cumsum([SaveDiffData(index),PreR]);endend%-------------------预测趋势并返回结果----------------mp1=polyfit([1:size(TrendData',2)],TrendData',1);xt=[];forj=1:stepxt=[xt,size(TrendData',2)+j];endTrendResult=polyval(mp1,xt);PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2));tempx=[TrendData',TrendResult]+PreR;%tempx为预测结果plot(tempx,'r-.');holdonplot(Data,'b');legend('ARMA拟合时序曲线','实际时序风速');save('resultarma.mat');