mk突变点检测_科学网—从网上找的M-K突变检验的程序 - 张乐乐的博文
%从matlab论坛上找的MK突变检验的程序,这个程序运行的结果跟我自己编写程序运行出来的结果一样,但是跟魏凤英老师书上的例子出图结果不一样
A=xlsread('test-mk.xlsx');
x=A(:,1);%时间序列
y=A(:,2);%径流数据列
N=length(y);
n=length(y);
% 正序列计算---------------------------------
% 定义累计量序列Sk,长度=y,初始值=0
Sk=zeros(size(y));
% 定义统计量UFk,长度=y,初始值=0
UFk=zeros(size(y));
% 定义Sk序列元素s
s = 0;
% i从2开始,因为根据统计量UFk公式,i=1时,Sk(1)、E(1)、Var(1)均为0
% 此时UFk无意义,因此公式中,令UFk(1)=0
for 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,初始值=0
y2=zeros(size(y));
% 定义逆序累计量序列Sk2,长度=y,初始值=0
Sk2=zeros(size(y));
% 定义逆序统计量UBk,长度=y,初始值=0
UBk=zeros(size(y));
% s归0
s=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)=0
for 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文件:f:test2.xls
% 目标表单:Sheet1
% 目标区域:UFk从A1开始,UBk2从B1开始
%xlswrite('f:test2.xls',UFk,'Sheet1','A1');
%xlswrite('f:test2.xls',UBk2,'Sheet1','B1');
figure(3)%画图
plot(x,UFk,'r-','linewidth',1.5);
hold on
plot(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 on
hold on
plot(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);
转载本文请联系原作者获取授权,同时请注明本文来自张乐乐科学网博客。
链接地址:http://wap.sciencenet.cn/blog-1103122-843162.html
上一篇:书中关于MK突变检验的问题
下一篇:给小魏师妹写的数据处理程序
mk突变点检测_科学网—从网上找的M-K突变检验的程序 - 张乐乐的博文相关推荐
- matlab 二维矩形函数,科学网—利用MATLAB对非矩形域实现二维插值 - 张乐乐的博文...
>> load('x1.mat'); >> load('y1.mat') >> load('T.mat'); >> load('Lon.mat'); & ...
- matlab 二维插值 验证,科学网-利用MATLAB对非矩形域实现二维插值-张乐乐的博文...
>> load('x1.mat'); >> load('y1.mat') >> load('T.mat'); >> load('Lon.mat'); & ...
- java swing 左上角图标_科学网—Matlab: 学习GUI(修改窗口左上角图标而不warning) - 刘磊的博文...
网上常用的方法: if ~isdeployed newIcon=javax.swing.ImageIcon('.piciap.jpg'); else newIcon=javax.swing.Image ...
- endnote文件enl突然没了_科学网—实际操作中的Endnote库文件损坏修复方法 - 尹卓忻的博文...
Endnote是保存文件的神器,将文献的详细信息输入标签之后,插入文献只用点一下.不过就算是神器也有掉链子的时候,有时内力不够,刚打开就跳出以下界面: 按对话框的信息,问题是可以通过重启恢复 , ...
- python 面板数据分析_科学网—Python中的结构化数据分析利器-Pandas简介 - 郑俊娟的博文...
此文转载于XXXXXX处... Pandas是python的一个数据分析包,最初由AQR Capital Management于2008年4月开发,并于2009年底开源出来,目前由专注于Python数 ...
- python序列_科学网—Python:序列(字符串、列表、元组)和序列函数 - 刘洋洋的博文...
Python中的序列,包括字符串(String).列表(List).元组(Tuple). 序列的索引 通过索引(index)访问及获得的序列的一个或多个元素,也叫切片. 正序: 0 到 N-1 倒序: ...
- 探测器反向偏压_科学网—《涨知识啦22》---MSM型光电探测器 - 寇建权的博文
此前,小赛给大家简单普及了金属与半导体之间的两种接触类型:欧姆接触与肖特基接触,二者也凭借各自的优势被研究人员充分应用.本周小赛给大家主要介绍的是基于肖特基接触类型的MSM型光电探测器的基本原理. 众 ...
- matlab mic系数_科学网—最大信息系数 (Maximal Information Coefficient, MIC)详解(1) - 彭勇的博文...
最大信息系数 (Maximal Information Coefficient, MIC)详解(1) 四年前看过的一篇论文,当时还在组会上报告过,很确信当时把它弄懂了,由于当时是用机器学习的方法来做预 ...
- python读取tiff影像_科学网—利用python GDAL库读写geotiff格式的遥感影像方法 - 张伟的博文...
(1)利用python GDAL库读写geotiff格式的遥感影像方法,具有很好的参考价值,不错! from osgeo import gdal import numpy as np def read ...
最新文章
- lvs+keepalived+nginx+tomcat高可用高性能集群部署
- 作用于HTML元素的Vue.js指令
- 不要再被Python洗脑了,来看看这个吧......
- LR回归原理和损失函数的推导
- 李子奈计量经济学笔记和课后习题答案
- openwrt 添加usb网卡_如何在双网卡的蜗牛星际上打一个openwrt软路由和NAS一体机
- python爬取中国大学排名_Python爬取中国大学排行榜
- QT中使用以管理员权限启动一个进程
- 感恩节特辑丨这个世界有无限可能
- 创建通用 macOS 二进制文件
- 如何删除计算机新用户,如何将电脑里的账户信息彻底删除
- 输入多个数,中间用空格隔开
- Mybatis关联查询的两种方式
- oracle怎么算时间,Oracle时间计算
- 山东大学软件工程应用与实践——WeaselTSF(一)
- 2017常见android面试题
- 笑抽了~~关于程序员的爆笑gif图片
- python爬取智联招聘网_python爬取智联招聘工作岗位信息
- 日化美妆难突围,看爱码物联如何冲破传统营销壁垒
- 岭南(含广东广西海南)地形及DEM下载
热门文章
- lay和lied_lay和lie的区别
- Android打包apk实现原理与流程(雷惊风)
- 卖计算机配件的二手平台,电脑哪些配件适合买二手,哪些最好入新?
- phpadmin 导入数据
- 列联表相关测量--c相关系数
- 山东高中学业水平考试时间2020计算机,2020年山东省高中学业水平等级考试报名时间及科目...
- 解决“error: failed to push some refs to ‘git@gitee.com:username/repo.git‘“
- html5怎么做相册影集,手机怎么做相册影集
- Node.js 第一天
- flask项目中出现Error: While importing ‘manager‘, an ImportError was raised.