LoginSignup
2
1

More than 1 year has passed since last update.

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

Last updated at Posted 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;
}

参考文献

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

関連論文

  1. ROBERT M. HIRSCH AND JAMES R. SLACK(1984), A Nonparametric Trend Test for SeasonalData With Serial Dependence, WATER RESOURCES RESEARCH, VOL. 20, NO. 6, PAGES 727-732
2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
1