利用R语言进行MK检验以检测序列变化趋势

关于MK检验，目前因为搜索不到相关的R语言程序包，所以本人按照魏凤英老师在其著作《现代气候统计诊断与预测技术》中的MK检验方法，利用R语言写下如下代码，希望能够帮到大家，相互切磋，相互交流。

本文所用数据也为魏老师书中1900-1990年上海年平均气温序列

注：本文只是为大家提供思路，文中所画图形与魏凤英老师书中对应图像有所差异，但我觉得的这个程序没有错，本人才疏学浅，如果看到这篇文章的您觉得我程序中哪里有错误，望您留言，吾不胜感激，多谢

#####设置文件路径######
setwd("D:\\paper\\data\\")
#####原始数据#####
library(openxlsx)
data<-read.xlsx("shanghai_testdata.xlsx",colNames = FALSE)
data<-as.matrix(data)
#####对数据求倒序#####
data_dx<-rev(data)
data_dx<-matrix(data_dx)
n<-length(data)
#######MK检验#######
Q<-matrix(0,1,n)
UF<-matrix(0,1,n)
UB<-matrix(0,1,n)
h<-matrix(0,n,n)
for (i in 1:n) {
  for (j in 1:n) {
    if(data[i,1]>data[j,1]){
      h[i,j]<-1
    }
    else
      h[i,j]<-0
  }
  Q[1,i]<-sum(h[lower.tri(h)])
  UF[1,i]<-(Q[1,i]-(i*(i-1)/4))/sqrt((i*(i-1)*(2*i+5))/72)
}
#####计算UB#####
h<-matrix(0,n,n)
for (i in 1:n) {
  for (j in 1:n) {
    if(data_dx[i,1]>data_dx[j,1]){
      h[i,j]<-1
    }
    else
      h[i,j]<-0
  }
  Q[1,i]<-sum(h[lower.tri(h)])
  UB[1,i]<-(Q[1,i]-(i*(i-1)/4))/sqrt((i*(i-1)*(2*i+5))/72)
}
#####绘图#####
UF[1,1]<-0
UB[1,1]<-0
plot(x=1900:1990,y=UF,ylim = c(-4,8),type = "l",ylab = "",xaxt="n")
lines(x=1900:1990,y=-rev(UB),type = "l",lty=2)
lines(x=1900:1990,y=rep(1.96,91),type = "l")
lines(x=1900:1990,y=rep(-1.96,91),type = "l")
lines(x=1900:1990,y=rep(0,91),type = "l")
title("MK检验_上海")
axis(1,1900:1990,1900:1990,las=1)
#####输出突变年份#####
year_mk<-1900:1990
year_point<-year_mk[which((as.numeric(UF)-(-rev(UB)))>0)[1]-1]
print(year_point)

输出突变年份结果如下：