Graviness Blog

JUGEMテーマ：インターネット

JUGEMテーマ：コンピュータ

JUGEMテーマ：経済全般

JUGEMテーマ：面白数学・面白科学

JUGEMテーマ：学問・学校

JUGEMテーマ：学問・学校

JUGEMテーマ：コンピュータ

スミルノフ−グラブス検定を用いる外れ値除去のPython（NumPy/SciPy使用）の実装です。

下記smirnov_grubbs関数にサンプルデータと有意水準α（両側100α％）を入力すると、外れ値を除去したサンプルデータと、抽出した外れ値データを出力します。

サンプルデータのうち平均値からの偏差が最も大きいデータを、有意水準（両側）でスミルノフ−グラブス検定し、帰無仮説を棄却した場合に当該データを外れ値とします。当該データを除いたサンプルデータに対し、外れ値がなくなるまでこの操作を繰り返します。

ref. スミルノフ−グラブス検定＠wikipedia

# coding: utf-8
#
# smirnov_grubbs(data, alpha)
#   Generate outlier-removed data with Smirnov-Grubbs.
#
#   Parameters:
#     data: array_like
#       Numeric data values.
#     alpha: float
#       Significance level. (two-sided)
#
#   Returns:
#     outlier-removed-data
#     outlier-data
#
#   e.g.
#     data = [1, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 1000]
#     alpha = 0.01
#     smirnov_grubbs(data, alpha)
#     ==> ([100, 101, 102, 103, 104, 105, 106, 107, 108, 109], [1000, 1])

import numpy as np
import scipy.stats as stats

def smirnov_grubbs(data, alpha):
    x, o = list(data), []
    while True:
        n = len(x)
        t = stats.t.isf(q=(alpha / n) / 2, df=n - 2)
        tau = (n - 1) * t / np.sqrt(n * (n - 2) + n * t * t)
        i_min, i_max = np.argmin(x), np.argmax(x)
        myu, std = np.mean(x), np.std(x, ddof=1)
        i_far = i_max if np.abs(x[i_max] - myu) > np.abs(x[i_min] - myu) else i_min
        tau_far = np.abs((x[i_far] - myu) / std)
        if tau_far < tau: break
        o.append(x.pop(i_far))
    return (np.array(x), np.array(o))

JUGEMテーマ：コンピュータ

Brunner-Munzel検定のPython（NumPy/SciPy使用）の実装です。

下記bmtest関数に2群のデータを入力すると、検定統計量W、p値（両側）、自由度、および帰無仮説で示される確率の推定値を出力します。

片側検定にする場合（どちらにずれているかが重要な場合）、出力されたp値を0.5倍して、Wの値が正負で判断します。

注意: 少しもかぶりもしない2群を入力すると、分散が0となって、検定統計量Wや自由度の計算中に Runtimewarning; divide by zero が出ます。

ref. Brunner-Munzel検定＠奥村 晴彦さん

# coding: utf-8
#
# bmtest(x1, x2)
#   Calculate Brunner-Munzel-test scores.
#
#   Parameters:
#     x1, x2: array_like
#       Numeric data values from sample 1, 2.
#
#   Returns:
#     w:
#       Calculated test statistic.
#     p_value:
#       Two-tailed p-value of test.
#     dof:
#       Degree of freedom.
#     p:
#       "P(x1 < x2) + 0.5 P(x1 = x2)" estimates.
#
#   References:
#     * https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

import numpy as np
import scipy.stats as stats

def bmtest(x1, x2):
    n1, n2 = len(x1), len(x2)
    R = stats.rankdata(list(x1) + list(x2))
    R1, R2 = R[:n1], R[n1:]
    r1_mean, r2_mean = np.mean(R1), np.mean(R2)
    Ri1, Ri2 = stats.rankdata(x1), stats.rankdata(x2)
    var1 = np.var([r - ri for r, ri in zip(R1, Ri1)], ddof=1)
    var2 = np.var([r - ri for r, ri in zip(R2, Ri2)], ddof=1)
    w = ((n1 * n2) * (r2_mean - r1_mean)) / ((n1 + n2) * np.sqrt(n1 * var1 + n2 * var2))
    dof = (n1 * var1 + n2 * var2) ** 2 / ((n1 * var1) ** 2 / (n1 - 1) + (n2 * var2) ** 2 / (n2 - 1))
    c = stats.t.cdf(abs(w), dof) if not np.isinf(w) else 0.0
    p_value = min(c, 1.0 - c) * 2.0
    p = (r2_mean - r1_mean) / (n1 + n2) + 0.5
    return (w, p_value, dof, p)

JUGEMテーマ：コンピュータ