matlab程序.docx
《matlab程序.docx》由会员分享，可在线阅读，更多相关《matlab程序.docx（14页珍藏版）》请在冰豆网上搜索。
matlab程序
clear
clc
fid=fopen（'CO2.txt','r'）;
tline=fgetl（fid）;
co=fscanf（fid,'%f',[13,57]）;
Fa=2;
n=1;
fori=1:
1:
53%把矩阵改写为时间序列
ifi==1
forj=4:
1:
13
tco（n）=co（j,i）;
t（n）=co（1,i）+j./12;
n=n+1;
end
else
forj=2:
1:
13
tco（n）=co（j,i）;
t（n）=co（1,i）+j./12;
n=n+1;
end
end
end
%%%
%%趋势项处理
X=ones（n-1,2）;
X（:
2）=t';
b=regress（tco',X）;
p2=polyfit（t,tco,2）;
fori=1:
1:
n-1
q1（i）=b
（1）+b
（2）*t（i）;%线性拟合
q2（i）=p2
（1）*t（i）.^2+p2
（2）*t（i）+p2（3）;%二次拟合
end
figure
（1）
plot（t,tco）
title（'CO_{2}浓度的时间变化曲线','fontsize',15）
xlabel（'Time','fontsize',15）
ylabel（'CO_{2}浓度/ppm','fontsize',15）
holdon
%plot（t,q1,'linestyle','--','color',[100]）
plot（t,q2,'linestyle','-','color',[100]）
legend（'实测曲线','趋势线','location','northwest'）
legend（'boxoff'）
%%%
%%周期项处理
%%这里用谐波分析，假设T=634，这样处理有问题但结果没有受到太大影响
fori=1:
1:
n-1
c（i）=tco（i）-q2（i）;%剔除趋势项
end
cc=smooth（c,12）;%12点平滑
a0=mean（c）;
w=2*pi/（n-1）;
fori=1:
1:
floor（（n-1）/2）
forj=1:
1:
n-1
ccc（i,j）=cos（i*w*j）;
sss（i,j）=sin（i*w*j）;
end
a（i）=2*mean（c.*ccc（i,:
））;
b（i）=2*mean（c.*sss（i,:
））;
s（i）=（a（i）.^2+b（i）.^2）/2;%计算方差贡献
end
i=1:
1:
floor（（n-1）/2）;
figure
（2）
plot（i,s）
title（'谐波谱图','fontsize',15）
xlabel（'波数','fontsize',15）
ylabel（'方差贡献','fontsize',15）
[news,I]=sort（s,'descend'）;%方差贡献从小到大排列
F=0;
kk=1;
whilekk<=length（s）%找显著波数
F=s（kk）*（n-4）/2/（var（c）-s（kk））;
ifF>Fa
kk=kk+1;
else
break
end
end
fori=1:
1:
n-1
p（i）=a0;
forj=1:
1:
kk-1
p（i）=p（i）+a（I（j））*cos（I（j）*w*i）+b（I（j））*sin（I（j）*w*i）;%周期项拟合
ppppp=（n-1）/I（j）%把波数变为周期
end
end
figure（3）
plot（t,c）
xlabel（'Time','fontsize',15）
ylabel（'CO_{2}浓度/ppm','fontsize',15）
holdon
plot（t,cc,'linestyle','-','color',[000]）
plot（t,p,'linestyle','--','color',[100]）
legend（'剔除趋势项后的变化曲线','12点平滑后的曲线','拟合曲线','location','northeast'）
legend（'boxoff'）
%%%
%%随机项处理
fori=1:
1:
n-1
sj（i）=c（i）-p（i）;%剔除周期项
end
u=mean（sj）;
fori=1:
1:
n-1
sj（i）=sj（i）-u;%得到平稳距平序列
end
figure（4）
i=1:
1:
n-1;
plot（i,sj）
title（'平稳随机序列曲线','fontsize',15）
xlabel（'Time','fontsize',15）
ylabel（'CO_{2}浓度/ppm','fontsize',15）
forto=1:
1:
floor（n/4）
K（to）=0;
fori=1:
1:
n-1-to
K（to）=K（to）+（sj（i）-mean（sj））*（sj（i+to）-mean（sj））;
end
K（to）=K（to）/（n-1-to）;
ro（to）=K（to）/（var（sj））;%计算相关函数，to是时间差
end
fori=1:
1:
floor（n/4）
forj=1:
1:
floor（n/4）
ifi==j
xro（i,j）=1;
elseifixro（i,j）=ro（j-i）;
else
xro（i,j）=ro（i-j）;
end
end
yro（i）=ro（i）;
end
fork=1:
1:
floor（n/4）
fi=xro（1:
k,1:
k）\yro（1:
k）';
FPE（k）=0;
fori=1:
1:
k
FPE（k）=FPE（k）+fi（i）*ro（i）;
end
FPE（k）=（1-FPE（k））*（n+k）/（n-k-2）;
end
[newFPE,I2]=sort（FPE,'descend'）;
zyp=I2（length（FPE））;%定阶，可能有问题
fi=xro（1:
zyp,1:
zyp）\yro（1:
zyp）';%计算自回归模型的系数
%%%
%%预报
nn=n;
fori=54:
1:
57%把2011-2014的数据改写为序列形式
ifi==57
forj=2:
1:
12
tco（nn）=co（j,i）;
t（nn）=co（1,i）+j./12;
nn=nn+1;
end
else
forj=2:
1:
13
tco（nn）=co（j,i）;
t（nn）=co（1,i）+j./12;
nn=nn+1;
end
end
end
fori=n:
1:
nn-1%计算随机项的拟合值
sj（i）=0;
forj=i-zyp:
1:
i-1
sj（i）=sj（i）+fi（i-j）*sj（j）;
end
sj（i）=sj（i）+u;
end
fori=n:
1:
nn-1
p（i）=a0;
forj=1:
1:
kk-1
p（i）=p（i）+a（I（j））*cos（I（j）*w*i）+b（I（j））*sin（I（j）*w*i）;
end
fco（i）=p2
（1）*t（i）.^2+p2
（2）*t（i）+p2（3）+p（i）+sj（i）;%预报结果
end
figure（5）
scatter（tco（n:
nn-1）,fco（n:
nn-1））
title（'实测值与拟合值散点图','fontsize',15）
xlabel（'实测值/ppm','fontsize',15）
ylabel（'拟合值/ppm','fontsize',15）
xx=tco（n:
nn-1）;
yy=fco（n:
nn-1）;
r=corrcoef（xx',yy'）%计算相关系数
figure（6）
plot（t（n:
nn-1）,tco（n:
nn-1））
holdon
plot（t（n:
nn-1）,fco（n:
nn-1）,'linestyle','--','color',[100]）
title（'实测曲线与拟合曲线','fontsize',15）
xlabel（'Time','fontsize',15）
ylabel（'CO_{2}浓度/ppm','fontsize',15）
legend（'实测曲线','拟合曲线','location','northwest'）
legend（'boxoff'）
花粉
fid=fopen（'D:
\gongju\matlab\mfile\dibiaoshuju.txt','r'）;
data=fscanf（fid,'%f',[10,74]）;
data=data';
fclose（fid）;
fid=fopen（'D:
\gongju\matlab\mfile\ttt.txt','r'）;
data3=fscanf（fid,'%f',[5,148]）;
fclose（fid）;
time=data3（1,:
）';
x=data3（2:
5,:
）';
x=[ones（148,1）x];
b=[10.1093-0.011847-0.079571-0.035224-0.015373]';
y=x*b;
b1=[1456.25790.370572.44331.15970.45866]';
y1=x*b1;
plot（time,y）;
title（'温度变化'）;
saveas（gcf,['D:
\gongju\matlab\mfile\tu\temperature.jpg']）;
plot（time,y1）;
title（'降水量变化'）;
saveas（gcf,['D:
\gongju\matlab\mfile\tu\rainfall.jpg']）;
fid=fopen（'D:
\gongju\matlab\mfile\height.txt','r'）;
height=fscanf（fid,'%f'）;
fclose（fid）;
fid=fopen（'D:
\gongju\matlab\mfile\tp.txt','r'）;
data2=fscanf（fid,'%f',[3,43]）;
fclose（fid）;
fori=1:
74
forj=1:
42
if（（height（i）>data2（1,j））&&（height（i）<=data2（1,j+1）））||（（j==42）&&（height（i）>data2（1,43）））||（（j==1）&&（height（i）<=data2（1,1）））
t（i）=data2（2,j）+（height（i）-data2（1,j））/（data2（1,j+1）-data2（1,j））*（data2（2,j+1）-data2（2,j））;
p（i）=data2（3,j）+（height（i）-data2（1,j））/（data2（1,j+1）-data2（1,j））*（data2（3,j+1）-data2（3,j））;
end
end
end
t=t';
p=p';
逐步回归
functionstepregress（x,y,F）
%x=zscore（x,1）;
%y=zscore（y,1）;
r=corrcoef（[x,y]）;
l=0;
L=0;
[n,m]=size（x）;
k=ones（m）;
q=1;
while（q==1）
q=0;
fori=1:
m
v（i）=r（i,m+1）^2/r（i,i）;
end
max=1;
min=1;
fori=1:
m
if（（max==1）&&（k（i）==1）&&（k
（1）==0））||（（v（i）>v（max））&&（k（i）==1））
max=i;
end
if（（min==1）&&（k（i）==0）&&（k
（1）==1））||（（v（i）min=i;
end
end
if（l<3）&&（L+1<=m）
F1=v（max）/（（r（m+1,m+1）-v（max））/（n-l-2））;
if（F1>F）
disp（['引入第',num2str（max）,'个变量']）;
k（max）=0;
L=L+1;
l=l+1;
r=matdel（max,m+1,r）;
q=1;
end
else
F2=v（min）/（r（m+1,m+1）/（n-l-1））;
if（（F2disp（['剔除第',num2str（min）,'个变量']）;
k（min）=1;
L=L-1;
l=l+1;
r=matdel（min,m+1,r）;
q=1;
else
F1=v（max）/（（r（m+1,m+1）-v（max））/（n-l-2））;
if（F1>F）
disp（['引入第',num2str（max）,'个变量']）;
k（max）=0;
L=L+1;
l=l+1;
r=matdel（max,m+1,r）;
q=1;
end
end
end
end
disp（'没有可剔除或引入的变量，逐步回归结束'）;
a=zeros（L）;
j=1;
fori=1:
m
if（k（i）==0）
a（j）=i;
j=j+1;
end;
end;
xx=x（:
a
（1））;
fori=2:
L
xx=[xxx（:
a（i））];
end;
xx=[ones（n,1）xx];
b=regress（y,xx）;%回归系数
R=sqrt（1-r（m+1,m+1））;%复相关系数
yyy=xx*b;%y的估计值
ymean=mean（y）;%y平均值
Q=（y-yyy）'*（y-yyy）;%剩余平方和
U=（yyy-ymean）'*（yyy-ymean）;%回归平方和
rs=Q/（n-L-1）;%剩余方差
f=U/L/（Q/（n-L-1））;%F统计量
fid=fopen（'result','w'）;
ss=['引入第',num2str（a
（1））];
fori=2:
L
ss=[ss,',',num2str（a（i））];
end
ss=[ss,'个变量'];
ss1=['y=',num2str（b
（1））,'+（',num2str（b
（2））,'x',num2str（a
（1））,'）'];
fori=2:
L
ss1=[ss1,'+（',num2str（b（i+1））,'x',num2str（a（i））,'）'];
end;
ss2=['复相关系数=',num2str（R）];
ss3=['剩余方差=',num2str（rs）];
ss4=['F统计量=',num2str（f）];
ss5=['剩余平方和=',num2str（Q）];
fprintf（fid,'%s\n',ss）;
fprintf（fid,'%s\n',ss1）;
fprintf（fid,'%s\n',ss2）;
fprintf（fid,'%s\n',ss3）;
fprintf（fid,'%s\n',ss4）;
fprintf（fid,'%s',ss5）;
fclose（fid）;
end
clear;
fid=fopen（'D:
\gongju\matlab\mfile\dicengshuju.txt','r'）;
data=fscanf（fid,'%f',[27,148]）;
t=data（1,:
）;
fork=1:
26
plot（t,data（k+1,:
））;
switchk
case1
title（'Pinus'）;
case2
title（'Tsuga'）;
case3
title（'Picea'）;
case4
title（'Betula'）;
case5
title（'Alnus'）;
case6
title（'Carpinus'）;
case7
title（'Corylus'）;
case8
title（'Quercus（落）'
end
saveas（gcf,['D:
\gongju\matlab\mfile\tu\tu',num2str（k）,'.jpg']）;
end;
升级会员 