爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 7281|回复: 13

[程序设计] m-k检验matlab程序出错,求大神指点!

[复制链接]

新浪微博达人勋

发表于 2014-5-12 19:54:57 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
% Mann-Kendall突变检测 % 数据序列y% 结果序列UFk,UBk2%--------------------------------------------%读取excel中的数据,赋给矩阵y%获取y的样本数%A为时间和气温数据列A=xlswrite('d:\data.xls')x=A(:,1);%时间序列y=A(:,2);%气温数据列N=length(y);n=length(y);% 正序列计算---------------------------------% 定义累计量序列Sk,长度=y,初始值=0Sk=zeros(size(y));% 定义统计量UFk,长度=y,初始值=0UFk=zeros(size(y));% 定义Sk序列元素ss = 0;% i从2开始,因为根据统计量UFk公式,i=1时,Sk(1)、E(1)、Var(1)均为0% 此时UFk无意义,因此公式中,令UFk(1)=0for i=2:n   for j=1:i         if y(i)>y(j)           s=s+1;         else           s=s+0;         end;   end;   Sk(i)=s;   E=i*(i-1)/4; % Sk(i)的均值  Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差  UFk(i)=(Sk(i)-E)/sqrt(Var);end;% ------------------------------正序列计算end% 逆序列计算---------------------------------% 构造逆序列y2,长度=y,初始值=0y2=zeros(size(y));% 定义逆序累计量序列Sk2,长度=y,初始值=0Sk2=zeros(size(y));% 定义逆序统计量UBk,长度=y,初始值=0UBk=zeros(size(y));% s归0s=0;% 按时间序列逆转样本y% 也可以使用y2=flipud(y);或者y2=flipdim(y,1);for i=1:n    y2(i)=y(n-i+1);end;% i从2开始,因为根据统计量UBk公式,i=1时,Sk2(1)、E(1)、Var(1)均为0% 此时UBk无意义,因此公式中,令UBk(1)=0for i=2:n   for j=1:i         if y2(i)>y2(j)           s=s+1;         else           s=s+0;         end;   end;   Sk2(i)=s;   E=i*(i-1)/4; % Sk2(i)的均值  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差% 由于对逆序序列的累计量Sk2的构建中,依然用的是累加法,即后者大于前者时s加1,% 则s的大小表征了一种上升的趋势的大小,而序列逆序以后,应当表现出与原序列相反% 的趋势表现,因此,用累加法统计Sk2序列,统计量公式(S(i)-E(i))/sqrt(Var(i))% 也不应改变,但统计量UBk应取相反数以表征正确的逆序序列的趋势  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);end;% ------------------------------逆序列计算end% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量% 与UFk做图寻找突变点时,2条曲线应具有同样的时间轴,因此% 再按时间序列逆转结果统计量UBk,得到时间正序的UBk2,做图用UBk2=zeros(size(y));% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);for i=1:n   UBk2(i)=UBk(n-i+1);end;% 做突变检测图时,使用UFk和UBk2% 写入目标xls文件:d:\test2.xls% 目标表单:Sheet1% 目标区域:UFk从A1开始,UBk2从B1开始B=xlswrite('d:\test2.xls',UFk,'Sheet1','A1');B=xlswrite('d:\test2.xls',UBk2,'Sheet1','B1');figure(3)%画图plot(x,UFk,'r-','linewidth',1.5);hold onplot(x,UBk2,'b-.','linewidth',1.5);plot(x,1.96*ones(N,1),':','linewidth',1);axis([min(x),max(x),-5,5]);legend('UF统计量','UB统计量','0.05显著水平');xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);%grid onhold onplot(x,0*ones(N,1),'-.','linewidth',1);plot(x,1.96*ones(N,1),':','linewidth',1);plot(x,-1.96*ones(N,1),':','linewidth',1)

出错提示为Error using ==> xlswriteNot enough input arguments.


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2014-5-12 19:56:31 | 显示全部楼层
大家还是不要看这个帖子了,我自己看了都受不了,不知道为什么会变成这个样子,好好的一个程序,唉!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-5-12 20:17:01 | 显示全部楼层
那是因为你没有排版好,直接摘过来的吧?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2014-5-12 21:01:27 | 显示全部楼层
kongfeng0824 发表于 2014-5-12 20:17
那是因为你没有排版好,直接摘过来的吧?

从txt粘过来的,不过问题已经解决了!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-5-13 08:55:25 | 显示全部楼层
DGT 发表于 2014-5-12 21:01
从txt粘过来的,不过问题已经解决了!

弱弱的问一句,您这个程序没有错误,可以直接用吗?我也要用M-K检验,我是新手
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-5-13 08:56:45 | 显示全部楼层
可以直接用么?我也要用这个检验,我是新手,不太懂
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2014-5-14 17:29:15 | 显示全部楼层
小鸽子1991125 发表于 2014-5-13 08:56
可以直接用么?我也要用这个检验,我是新手,不太懂

这个程序有个错,改对了应该可以直接出图,但我不知道错在哪,后来用别的方法做的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-7-25 10:50:15 | 显示全部楼层
也好想直接用啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-7-26 20:18:36 | 显示全部楼层
国标《农业干旱预警等级》(征求意见稿)附录B有标准
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-8-19 21:43:09 | 显示全部楼层
重新排版吧
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表