时间序列分析及MK检验的原理

时间序列分析是研究随时间变化而变化的现象和规律的一种分析方法。其包括分析序列的趋势、周期性、季节性、随机性等成分，并利用这些成分对未来进行预测。

MK检验（Mann-Kendall test）是一种常用的非参数检验方法，用于检验时间序列数据中是否存在趋势。其基本原理是比较样本中每个数据点与其前面的数据点之间的大小关系，通过计算得到一个统计量，判断其是否显著不同于零。

MK检验的步骤如下：

1.计算每个数据点与其前面的数据点之间的差值，得到一个序列。

2.计算序列中每个差值为正数或负数的排列组合数。

3.计算序列中所有差值的秩和，得到统计量。

4.对统计量进行正态近似或查表找出其对应的显著性水平。

5.根据显著性水平判断是否存在趋势。

需要注意的是，MK检验假设时间序列中的数据是独立同分布的，且不存在季节性和周期性等特征。如果时间序列中存在这些特征，需要进行相应的预处理或采用其他方法进行趋势分析。

相关问题

多年NDVI数据时间序列数据的MK检验matlab

### 使用Matlab对多年NDVI时间序列数据执行Mann-Kendall趋势检验

#### 准备工作
为了在Matlab中进行Mann-Kendall趋势检验，首先需要准备NDVI时间序列数据。假设这些数据已经整理成一个向量形式的时间序列`ndvi_series`。

#### 实现步骤概述
可以利用Matlab中的函数来完成这一任务。虽然Matlab官方库未直接提供Mann-Kendall测试功能，但是可以通过编写自定义函数或使用社区贡献的工具箱来进行此操作。下面给出一种基于已有资源的方法：

#### 自定义Mann-Kendall Test Function
如果希望不依赖额外包而实现基本的功能，则可以根据算法原理构建自己的函数。这里展示了一个简单的例子：

function [tau, pval] = mann_kendall(x)
    % 计算S统计量
    n = length(x);
    S = 0;
    for i=1:n-1
        for j=i+1:n
            if x(j)>x(i)
                S = S + 1;
            elseif x(j)<x(i)
                S = S - 1;
            end
        end
    end
    
    % 方差Var(S)计算
    ties = unique(x); % 找到所有的不同值
    g = numel(ties);
    VarS = (n*(n-1)*(2*n+5))/18; % 初始方差估计
    tiecorrection = 0;

    for k=1:g
        tpk = sum(x==ties(k)); % 当前值得重复次数
        if(tpk>1)
           tiecorrection = tiecorrection+(tpk*(tpk-1)*(2*tpk+5));
        end
    end
    VarS = VarS-tiecorrection/18;

    % 正态近似得到Z分数
    Z = sign(S)*((abs(S)-1)/sqrt(VarS));

    % Kendall's Tau 和 P-value 的计算
    tau = S/(n*(n-1)/2);
    pval = 2 * normcdf(-abs(Z), 0, 1);

end

上述代码实现了最基本的Mann-Kendall检验逻辑[^1]。通过调用该函数并传入NDVI时间序列作为输入参数即可获得Tau系数以及对应的P值。

#### 调用示例
假设有名为`ndvi_data.mat`的数据文件存储着多年的NDVI数值，在命令窗口中加载数据后可以直接调用上面编写的mann_kendall()函数:

load('ndvi_data.mat'); % 加载包含 NDVI 时间序列变量 ndvi_series 的 .mat 文件
[tau, pvalue] = mann_kendall(ndvi_series);
disp(['Kendall''s Tau: ', num2str(tau)]);
disp(['p-value: ', num2str(pvalue)]);

这段脚本会读取预先保存好的NDVI数据集，并对其进行Mann-Kendall趋势分析，最后打印出结果。

#### 结果解释
- `tau`: 表征两个随机选取样本之间顺序关系的概率差异程度；
- `pvalue`: 显示所观察到的趋势是否显著的小概率事件发生的可能性大小。通常情况下，当p<0.05时认为存在明显变化趋势。

mk趋势分析和mk突变检验

### Mann-Kendall 趋势分析方法

Mann-Kendall (MK) 趋势检验是一种非参数统计检验方法，广泛应用于检测时间序列数据中的单调趋势。该方法的优点在于不依赖于特定的数据分布形式，并且对异常值具有较强的鲁棒性[^3]。

#### 计算过程
1. **构建序列**：给定一组观测数据 \( X_1, X_2, \ldots, X_n \)，其中 n 是样本数量。
2. **计算S统计量**：对于每一对 i 和 j （i < j），如果 \( X_j > X_i \)，则记为 +1；反之，则记为 -1 或 0（相等情况下）。最终 S 统计量等于所有这些差值的总和。

\[
S = \sum_{k=1}^{n-1}\sum_{j=k+1}^n sgn(X_j-X_k)
\]

3. **方差 Var(S)** 的计算：

\[
var(s)=\frac{1}{18}[n(n-1)(2n+5)-\sum_{p=1}^q t_p(t_p-1)(2t_p+5)]
\]

这里 q 表示相同数值组的数量，\( t_p \) 则表示第 p 组内重复次数。

4. **标准化 Z 值**：当 n 大于或等于 10 时，可以近似认为 MK 检验服从标准正态分布 N(0,1)。此时可利用下述公式求得 z 分数：

如果 \( S>0 \):

\[
Z=\frac{S-1}{\sqrt{\text {var }(s)}}
\]

否则:

\[
Z=\frac{S+1}{\sqrt{\text {var }(s)}}
\]

5. **显著水平判断**：通过查表获得对应置信度下的临界值来决定原假设是否成立。通常采用双侧检验方式，在 α 显著性水平上拒绝零假设意味着存在上升/下降的趋势。

from pymannkendall import original_test

result = original_test([your_data_here])
print(f'Tau: {result.tau}, P-value: {result.p}')

### Mann-Kendall 突变检验方法

除了上述用于评估整体趋势外，还可以借助 Pettitt's test 来识别潜在的变化点位置——即所谓的“突变”。此技术同样基于秩次变换原理，旨在寻找最有可能发生结构性改变的时间节点。

具体实现如下所示：

1. 对原始序列按照升序排列并赋予相应排名；
2. 构造辅助函数 U(i), 它反映了前 i 项累积差异情况；
3. 寻找使得 |U(i)| 达到最大绝对值得那个时刻作为可能断点；
4. 结合 K-S 检验或其他手段进一步验证所得结论的有效性和可靠性。

import numpy as np
from scipy.stats import norm


def pettit_test(x):
    """Pettit’s change point detection."""
    
    n = len(x)
    k = range(1, n)

    r = stats.rankdata(x).astype(int)
    u = [abs(sum(r[:kk])-(kk)*(n+1)/2.) for kk in k]
    u_star = max(u)
    du = abs(np.diff(x))
    w = sum(du)+((min(du))**2)/(2*sum(du))

    p_value = 2*(1-norm.cdf(u_star/w*np.sqrt(2/n)))
    cp_loc = int(k[list(u).index(u_star)])

    return {'change_point_location': cp_loc,
            'u_statistic': u_star,
            'w_statistic': w,
            'p_value': p_value}

# Example usage with synthetic data containing a known changepoint at index 50.
np.random.seed(42)
synthetic_series = list(np.random.normal(loc=-1., scale=.5, size=(50))) \
                  +list(np.random.normal(loc=1., scale=.5, size=(50)))

pettit_result = pettit_test(synthetic_series)
print(pettil_result)