More than 1 year has passed since last update.

Mann-Kendall 検定を実装してみる

Last updated at 2022-10-16Posted at 2022-10-10

はじめに

時系列データにトレンド（傾向）があるかどうかを検定する方法の一つとして、Mann-Kendall検定が知られている。特徴としては、ノンパラメトリック法である。つまり、検定するサンプルがどんな分布に従っているかは仮定しない。その他、季節性検定バージョンもあるらしい。

検定方法

検定する $n > 1$ 個のサンプルを ${x_1,x_2,\ldots,x_n}$ とし、帰無仮説 $H_0$ と対立仮説 $H_1$ を以下とする。
$H_0$：$n$ 個のサンプル $x_1, x_2, \ldots, x_n$ が独立で同一の確率分布にしたがう。
$H_1$：$n$ 個のサンプル $x_1, x_2, \ldots, x_n$ が同一の確率分布にしたがわない。

以下、検定手順を示す。

STEP.1 統計量Sの計算

$$
S=\sum_{i=1}^{n-1}\sum_{j = i+1}^{n}\mathrm{sgn}(x_j-x_i)\in\mathbb{Z}
$$

ここで、

$$
\mathrm{sgn}(x)=
\begin{cases}
1 & x > 0 \\
0 & x = 0 \\
-1 & x < 0
\end{cases}.
$$

STEP.2 Tied Group に分ける

時系列データ $(x_i)$ の順序統計量 $(y_i)$ をとする。単純に昇順でソートすればよい。値が等しいものは同じ Tied Group に入る。$i$ 番目の Tied Group の要素数を $t_i$ とし、Tied Group の個数（組数）を $m$ とする。

具体的には、順序統計量が
$$
y_1 < y_2 < y_3 = y_4 = y_5 < y_6 = y_7 < y_8
$$
の場合は、5つの Tied Group に分かれる：
$$
y_1\ |\ y_2\ |\ y_3,\ y_4,\ y_5\ |\ y_6,\ y_7\ |\ y_8.
$$
このときは、以下のような値をとる：
$$
m = 5,\\
t_1 = t_2 = 1,\ t_3 = 3,\ t_4 = 2,\ t_5 = 1.
$$

STEP.3 統計量 S の分散の計算

$$
\sigma^2=\mathrm{Var}[S]=N_S-T_S
$$
ここで、
$$
N_S=\frac{n(n-1)(2n+5)}{18},
$$
$$
T_S=\sum_{i=1}^m\frac{t_i(t_i-1)(2t_i+5)}{18}.
$$

ちなみに $E[S]=0$ である。
統計量 S によって、$S > 0$ のとき増加トレンド、$S < 0$ のとき減少トレンドがあると判断できる。

STEP.4 統計量 S の基準化した標準統計量 Z の計算

$$
Z=\frac{S-\mathrm{sgn}(S)}{\sigma}=
\begin{cases}
\frac{S-1}{\sigma} & S > 0 \\
0 & S = 0 \\
\frac{S+1}{\sigma} & S < 0
\end{cases}
$$

統計量 $S$ はほぼ正規分布に従うので、標準化（っぽい変換？）した標準統計量 $Z$ はほぼ正規分布 $N(0,1)$ に従う。$E[S]=0$ なのでほとんど $\sigma$ で割っているだけである。

STEP.5 Kendall の τの計算（任意）

Kendall の $\tau$ は Kendall の順位相関係数ともいい、順位の相関の強さを表す。この係数は統計量 $S$ と密接な次の関係がある：
$$
\tau = \frac{S}{D}
$$
ここで、
$$
D=\sqrt{\frac{1}{2}n(n-1)-\frac{1}{2}\sum_{i=1}^mt_i(t_i-1)}\sqrt{\frac{1}{2}n(n-1)}.
$$

定義から $D \geq 0$ なので、$S$ と $\tau$ の符号は一致する。つまり、増減トレンドは $\tau$ でも判断でき、文献によっては $\tau$ を判定に採用しているものもある。

STEP.6 Sen's Slope の計算（任意）

しばしば、Mann-Kendall 検定でトレンドの有無を検定し、Sen's Slope でトレンドの大きさ（傾き、変化率）を推定する方法が用いられる。

$$
b = \mathrm{median}\left(\frac{x_j - x_i}{j - i}\middle|1\leq i<j\leq n\right)$$

中央値をとる集合の元は変化率を表しているので、$b$ は変化率の推定値を表していることが分かる。

STEP.7 p 値の計算

$p$ 値は、両側検定であれば $P(|X|>|Z|)$ なので、標準正規分布の累積分布関数 $\mathrm{cdf}$ を利用して、以下の通り計算できる。
$$
p = 2 (1-\mathrm{cdf}(|Z|)).
$$

$Z$ は負の場合もあるので絶対値をとる。両側検定なので、2倍している。

STEP.8 検定結果

有意水準を $\alpha$ とすると、$p < \alpha$ のとき帰無仮説 $H_0$ が棄却される。棄却されないとトレンドの増減について何も言えない。

統計量の符号

以下に統計量の符号の関係を示す。

$S$	$<-1$	$-1$	$0$	$1$	$>1$
$Z$	$-$	$0$	$0$	$0$	$+$
$\tau$	$-$	$-$	$0$	$+$	$+$

以上の関係から $S$ と $\tau$ の符号は一致する。$S$ と $Z$ の場合は、$S=\pm 1$ のとき、$Z=0$ となるので、符号は若干異なる。$Z$ で増減を判定していたものもあるが、どうなんだろう？

注意点

Mann-Kendall 検定には、以下の第2の過誤が報告されている：

サンプルの期間の長さが短いために、傾向変動が検出できない。
経年変化が小さいために、傾向を検出できない。

サンプルコード

MannKendallTester.cs

/// <summary>
/// 検定する
/// </summary>
/// <param name="samples">サンプル</param>
/// <returns>検定結果を返す。</returns>
public MannKendallTestResult Test(IEnumerable<double> samples)
{
    var n = samples.Count();
    var n1 = n - 1;
    var nn1 = n * n1;
    var nn1h = nn1 * 0.5;

    // STEP.1 統計量Sを計算する
    var S = 0.0;
    var samples1 = samples.Take(n1).Select((value, index) => (value, index));
    var q = new List<double>();
    foreach (var (value1, index1) in samples1)
    {
        var samples2 = samples.Skip(index1 + 1).Take(n - index1).Select((value2, index2) => (value2, index2 + index1 + 1));
        foreach (var (value2, index2) in samples2)
        {

            var k = Sign(value2 - value1, AlmostZero);
            S += k;
            q.Add((value2 - value1) / (index2 - index1));
        }
    }

    // STEP.2-3 分散を計算する
    var tiedGroupElementsCount = new List<int>();
    var orderedSamples = samples.OrderBy(s => s);
    var tmpSample = samples.First();
    var tmpCount = 0;
    foreach (var sample in orderedSamples.Skip(1))
    {
        if (Sign(sample - tmpSample, AlmostZero) == 0)
        {
            tmpCount++;
        }
        else
        {
            if (tmpCount > 0) { tiedGroupElementsCount.Add(tmpCount); }
            tmpSample = sample;
            tmpCount = 0;
        }
    }

    var N = nn1 * (2 * n + 5) / 18.0;
    var T = 0.0;
    var DT = 0.0;
    foreach (var t in tiedGroupElementsCount)
    {
        var tt = t * (t - 1);
        T += tt * (2 * t + 5);
        DT += tt;
    }
    T /= 18.0;
    DT *= 0.5;
    var V = N - T;

    // STEP.4 統計量Zを計算する
    var Z = 0.0;
    if (S > 0) { Z = (S - 1) / Math.Sqrt(V); }
    else if (S < 0) { Z = (S + 1) / Math.Sqrt(V); }

    // STEP.5 τを計算する
    var D = Math.Sqrt(nn1h - DT) * Math.Sqrt(nn1h);
    var tau = S / D;

    // STEP.6 スロープを計算する
    var slope = SummaryStatistics.Median(q);

    // STEP.7 p値を計算する
    var pValue = 2 * (1.0 - NormalDistribution.Cdf(Math.Abs(Z), 0.0, 1.0));
    var rejected = (pValue < SignificanceLevel);

    // STEP.8 傾向を求める
    var trend = TrendType.NoTrend;
    if (rejected)
    {
        if (S > 0) { trend = TrendType.Increasing; }
        else if (S < 0) { trend = TrendType.Decreasing; }
    }

    var result = new MannKendallTestResult()
    {
        PValue = pValue,
        Statistic = Z,
        Rejected = rejected,
        Trend = trend,
        Score = S,
        Tau = tau,
        VarianceOfScore = V,
        Slope = slope,
    };

    return result;
}

参考文献

西岡昌秋,寶馨(2003). Mann-Kendall 検定による水文時系列の傾向変動,Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 46 B
小林健一郎,宝馨,中北英一(2010). 全球気候モデル出力を用いた日本域の100 年確率日降水量の将来予測, 水工学論文集,第54巻
Thorsten Pohlert(2020). Non-Parametric Trend Tests and
Change-Point Detection
http://koshiba.sakura.ne.jp/math/stats/mannkendall.html
https://www.xlstat.com/ja/solutions/features/mann-kendall-trend-tests
https://vsp.pnnl.gov/help/vsample/design_trend_mann_kendall.htm
http://www7.civil.kyushu-u.ac.jp/kankyo/kaigan/2014/7.pdf
https://ja.wikipedia.org/wiki/ケンドールの順位相関係数
https://betashort-lab.com/データサイエンス/統計学/mann-kendallのトレンド検定/