分享个自己编的Matlab做全球格点的多元回归程序，外加导出dat数据

skylark

代码已经进行过检验了，没有问题！

另外，刚刚学Matlab一个月，自己编的程序还是有些复杂繁琐，大神们提提意见，新手大家一起共勉。

目的:做全球144x72格点的时间尺度为37*12=444的回归,每个格点的回归系数都不一样，但每个格点只有一个回归系数，这个回归系数是基于444个月份算出来的。注意在引用多个因子时要先确定因子之间的要基本独立，具体可以查matlab的regress函数对X的要求。

附上代码：

%Importdata

fid=fopen('e:/uh/ci.dat','rb','ieee-le');

na=fread(fid,inf,'float');

na=na(1:444); %取特定列数组中的一段

nb=load('e:/uh/n.txt'); %读取txt文件

pra=load('e:/uh/ph_pr_are.txt'); %grads对某点降水异常做的面积平均

prep=ncread('e:/uh/precip.mon.mean.nc','precip');% 读取nc数据

prep=prep(:,:,1:444); %提取特定时间段

prep=reshape(prep,10368,12,37);

a=mean(prep,3); % 对全球144x72个格点分别做1-12月的37年平均

%pp是气候态,37years

pp=cat(3,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a);  

prepa(:,:,:)=prep(:,:,:)-pp(:,:,:);% 全球降水异常

prepa1=reshape(prepa,144,72,444);

%根据corr选择相关最好的phase

nc=nb(:).*reshape(w(7:450),444,1);

%开始做回归%

X=[ones(size(nc)) nc na];

for i=1:144

for j=1:72

Y1=prepa1(i,j,1:444);

Y=reshape(Y1,444,1);

[b,bint,r,rint,stats] =regress(Y,X);

const(i,j)=b(1);

b_nc(i,j)=b(2);

b_na(i,j)=b(3);

end;

end;

%用回归系数算出重建场

for i=1:144

for j=1:72

for k=1:444

repra(i,j,k)=const(i,j)+b_nc(i,j)*nc(k)+b_na(i,j)*na(k);

end;

end;

end;

% 重建后的降水场

pp1=reshape(pp,144,72,444);

for i=1:144

for j=1:72

for k=1:444

repr(i,j,k)=repra(i,j,k)+pp1(i,j,k);

end;

end;

end;

%将数据输出dat文件格式

fid=fopen('e:/uh/repr1.dat','w');

for i=1:444

for j=1:72

for i=1:144

fwrite (fid,repra(i,j,k),'float32');

end;

end;

end;

fclose(fid);

二爷名声在外

思路清晰，程序还有很大的简化空间
在matlab里面能用矩阵计算就用矩阵计算，尽量少用循环

skylark

二爷名声在外 发表于 2017-3-29 17:06
思路清晰，程序还有很大的简化空间
在matlab里面能用矩阵计算就用矩阵计算，尽量少用循环

激动了好一会！版主回复我了！！确实很多地方没想到怎么用matrix解决，为了快速出结果就用循环了,刚刚入门，觉得Reshape把矩阵完美利用了！继续努力！

zihan6688

学习了，谢谢

fabula

挺不错的。

lgz

挺不错的。

二爷名声在外

skylark 发表于 2017-3-29 17:16
激动了好一会！版主回复我了！！确实很多地方没想到怎么用matrix解决，为了快速出结果就用循环了,刚刚入 ...

继续加油，看好你！

skylark

lgz 发表于 2017-3-29 20:18
挺不错的。

谢谢鼓励

skylark

zihan6688 发表于 2017-3-29 17:50
学习了，谢谢

共勉~

lancerv8

楼主加油！！！