Mann Kendall 突变检验及绘图
本帖最后由 zhouyuanshen 于 2014-5-28 13:48 编辑【作者】(必填)：青易
【问题描述】(必填)：进行非参数Mann Kendall突变检验,及绘图
【R语言代码】Mann_Kendall <- function(timeserial){#Mann Kendall 突变检验，传递参数
Mann_Kendall_sub <- function(timeserial){#需要进行两次秩的分析，因此在函数中嵌套了一个函数
r <- c()#分析的三个变量，具体含义可以参照魏凤英老师的书
s <- c()#秩序列。
U <- c()
for(i in 2:length(timeserial))#进行大小比较，从第二个开始与以前的数据进行大小比较
{r <- 0
for(j in 1:i)
{
if(timeserial>timeserial){r=r+1}#如果后面的数大于前面的数，则秩加1
}
s <- 0
for (ii in 2:i){
s <- s+r#秩序列。Sk是第i时刻数值大于ii时刻数值个数的累计数
}
U <- 0
U <- (s- (i*(i-1)/4))/sqrt(i*(i-1)*(2*i+5)/72)
}
r <- 0
s <- 0
U <- 0
LST <- list(r = r, s = s, U = U)
return (LST)
}
timeserial_rev <- rev(timeserial)#生成逆序列
y1 <- Mann_Kendall_sub(timeserial)#计算正序列
y2 <- Mann_Kendall_sub(timeserial_rev)#计算逆序列
y2$U <- -(rev(y2$U))#转换符号与顺序
LST <- list(UF=y1,UB=y2)
return(LST)
}
#这里是你要修改的地方
setwd("d:/")
od <- read.table("1.txt", header=T)
Variable <- c("Jan","Feb","Mar","Apr","May","Jun")
#修改上面代码
#可以自己定义，或者根据数据自动生成
rows <- length(Variable)
startyear <- as.numeric(od)
years <- od$Year
#如果要生成图片，请执行相应代码
tiff("filename.tif", width=14.6, height=16, units="cm", res=300, family = 'serif')
par(mfrow=c(rows,2),oma=c(3,0,0,0), mar=c(0,2,0,0),cex=0.7)
for(i in 1:length(Variable)){
name <- paste("od[ DISCUZ_CODE_0 ]quot;,Variable, sep="")
value <- eval(parse(text=name))
plot(value,type="l", ylab=Variable, cex.axis=0.6,xaxt="n",mgp=c(1,0.1,0),tck=-0.02)
if(i==length(Variable)){
axis(side=1, at=years, tck=-0.04, hadj=0.4, labels=years,mgp=c(1,0.4,0), cex.axis=1) # add x-axis to the last figure
axis(side=1, at=1:length(od$Year), tck=-0.01, hadj=0.4, labels=NA, cex.axis=1) # add month labels to the x-axis
mtext("年份",side=1,line=1.5)
}
d<-Mann_Kendall(value)#进行突变检验
yUF <- as.data.frame(d$UF)$U
yUB <- as.data.frame(d$UB)$U
plot(x=c(1:length(od$Year)),y=yUF, type="l", ylim=c(min(yUF,yUB,-1.96),max(yUF,yUB,1.96)),lwd=1, lty=5, ylab="", cex=0.5,xaxt="n",mgp=c(1,0.1,0),tck=-0.02)
points(x=c(1:length(od$Year)),y=yUB,type="l",lty=3,col=6,lwd=1)
abline(h=1.969,lty=4,lwd=0.5)# 1.969是a=0.05的显著性水平
abline(h=-1.96,lty=4,lwd=0.5)
abline(h=0,col="gray",lwd=0.5)
}
axis(side=1, at=years, tck=-0.04, hadj=0.4, labels=years,mgp=c(1,0.4,0), cex.axis=1) # add x-axis to the last figure
axis(side=1, at=1:length(od$Year), tck=-0.01, hadj=0.4, labels=NA, cex.axis=1) # add month labels to the x-axis
mtext("年份",side=1,line=1.5)
dev.off()【代码说明】：需改根据自身使用修改部分代码,已经标注了
【运行结果】：
【参考出处】：魏凤英的书
不错。。。欢迎楼主再接再厉 初学者，拜读下，呵呵 jonas 发表于 2014-6-14 21:52
不错。。。欢迎楼主再接再厉
谢谢鼓励！
stevanwei 发表于 2014-6-22 12:42
初学者，拜读下，呵呵
以后多交流！
楼主太给力了！就在找MK检验的代码，自己写了半天，最后用fortran搞定的。。。为毛R语言没有自己的包可以画图啊！？ 本帖最后由 zhuimengren 于 2015-11-11 17:06 编辑
非常好，赞一个 在实际运行的时候，数据会有一部分是缺失值，这样的话程序就无法运行了啊。这个在程序中如何处理呢？ 在实际运行的时候，数据会有一部分是缺失值，这样的话程序就无法运行了啊，总是提示Error in if (x > x) { : missing value where TRUE/FALSE needed
这个在程序中如何处理呢？ 大神啊，给你跪了！
页:
[1]
2