英招

牢骚太盛防肠断，风物长宜放眼量

0%

Mann-Kendall趋势突变检验(by R)

发表于 更新于 分类于 Valine：
通过R创建数学计算函数实现MannKendall趋势、突变检验
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#趋势检验函数
MannKendall_trend<- function(x) {
  n <- length(x)
  S <- 0
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      S <- S + sign(x[j] - x[i])
    }
  }
  varS <- (n * (n - 1) * (2 * n + 5)) / 18
  z <- (S - 1) / sqrt(varS)
  return(z)
}

#突变检验函数
Mann_Kendall_sudden <- function(timeserial){#Mann Kendall 突变检验，传递参数
  Mann_Kendall_sub <- function(timeserial){#需要进行两次秩的分析，因此在函数中嵌套了一个函数
    r <- c()#分析的三个变量
    s <- c()#秩序列。
    U <- c()
    
    for(i in 2:length(timeserial))#进行大小比较，从第二个开始与以前的数据进行大小比较
    {r[i] <- 0
    
    for(j in 1:i)
    {
      if(timeserial[i]>timeserial[j]){r[i] <- r[i]+1}#如果后面的数大于前面的数，则秩加1
    }
    
    
    s[i] <- 0
    for (ii in 2:i){
      s[i] <- s[i] + r[ii]#秩序列。Sk是第i时刻数值大于ii时刻数值个数的累计数
    }
    
    
    U[i] <- 0
    U[i] <- (s[i] - (i * (i - 1) / 4))/sqrt(i * (i - 1) * (2 * i + 5) / 72)
    
    }
    
    r[1] <- 0
    s[1] <- 0
    U[1] <- 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)  
}

#导入文件夹内的所有txt数据
txt_files <- list.files(path = "E:\\data\\txt_path", pattern = "\\.txt$", full.names = TRUE)  
txt_data_list <- lapply(txt_files, function(file) {  
  read.table(file, sep = "\t", header = TRUE)  
})

#趋势检验
for (i in 1:length(txt_files)) {
  print(txt_files[i])
  data1<-txt_data_list[[i]]
  colnames(data1)[4]<-"target"
  data<-data1$target
  data_del_na<-na.omit(data)
  
  print(MannKendall_trend(data_del_na))
}

#突变检验

for (i in 1:length(txt_files)) {
  i<-i+1
  print(txt_files[i])
  data1<-txt_data_list[[i]]
  colnames(data1)[4]<-"target"
  data<-data1$target
  cat(which(is.na(data), arr.ind=TRUE))
  data_del_na<-na.omit(data)
  
  result2 <- Mann_Kendall_sudden(data_del_na)  
  yUF <- as.data.frame(result2$UF[3])$U
  yUB <- as.data.frame(result2$UB[3])$U
  for (j in 1:(length(yUF)-1)) {
    if((yUF[j]-yUB[j])*(yUF[j+1]-yUB[j+1])<0){
      cat("突变点:",j)
    }
  }
}


#画图，结果保存在yUF、yUB，可导出用Origin画图不在R画
plot(x=c(1:length(yUF)),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(yUF)),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)