Skip to content

Instantly share code, notes, and snippets.

@c-bata
Created July 28, 2015 05:02

Revisions

  1. c-bata created this gist Jul 28, 2015.
    199 changes: 199 additions & 0 deletions SmirnovGrubbsTest.ipynb
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,199 @@
    {
    "cells": [
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "# スミルノフ・グラブス検定\n",
    "\n",
    "外れ値の検出手法として広く利用されている、スミルノフグラブス検定(Smirnov Grubbs Test)について整理。"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "## 標準偏差による検出\n",
    "\n",
    "外れ値検出というとまず思いつくのは、[外れ値 - Wiki](https://ja.wikipedia.org/wiki/%E5%A4%96%E3%82%8C%E5%80%A4) の一番上に載っているように、統計検定量を$\\frac{x - \\mu}{\\sigma}$とし、その値が有意点(2, 3等)より大きければ外れ値とみなす手法"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "## スミルノフグラブス検定による棄却\n",
    "\n",
    "スミルノフグラブス検定では、正規分布を仮定して、有意点を\n",
    "\n",
    "$\\frac{N-1}{\\sqrt{N}} \\times \\sqrt{\\frac{t^2}{N-2 + (t^2)}}$\n",
    "\n",
    "に設定する。ここで、一番大きな値(or 小さな値, or その両方)の統計検定量$\\frac{x - \\mu}{\\sigma}$ と有意点を比較し、棄却するかどうかを決める。"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "### t値の求め方\n",
    "\n",
    "有意点の式を見れば分かるようにt値が必要。\n",
    "t値を手計算で求めるなら [t分布表](http://www.koka.ac.jp/morigiwa/sjs/td.htm) とかを参考に求めるけど、今回はPythonのプログラムから求める。\n",
    "\n",
    "> ※ t値 を求める際には、片側検定と両側検定を混同しないように注意。\n",
    "理由がない限り両側検定を行う。\n",
    "\n",
    "例えば上のt分布表を見ると「自由度=8」の場合は1%水準(両側検定)では **2.3060** 。\n",
    "これを `scipy.stats` を用いて求める"
    ]
    },
    {
    "cell_type": "code",
    "execution_count": 16,
    "metadata": {
    "collapsed": false
    },
    "outputs": [],
    "source": [
    "from scipy import stats\n",
    "n = 8 # 自由度\n",
    "alpha = 0.05 # 危険率(この値に根拠はない。分野によるけど0.05ぐらいに設定することが多い。)"
    ]
    },
    {
    "cell_type": "code",
    "execution_count": 17,
    "metadata": {
    "collapsed": false
    },
    "outputs": [
    {
    "name": "stdout",
    "output_type": "stream",
    "text": [
    "1.85954803752\n"
    ]
    }
    ],
    "source": [
    "t = stats.t.isf(alpha, n)\n",
    "print(t)"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "こんな感じで, `scipy.stats.t.isf` に 危険率とデータ数(自由度)を与えれば t値が求まる。\n",
    "但し **2.3060** にならないといけないのに、値が違う。\n",
    "∵ 今回は両側検定だから。両側検定でやる際には、危険率を2で割る必要がある。\n",
    "理由については、下の図が分かりやすい。\n",
    "\n",
    "![片側検定と両側検定](https://dl.dropboxusercontent.com/spa/pxtopj46wcrqpk1/_vs57k-4.png)\n",
    "\n",
    "出展: http://www3.ocn.ne.jp/~stat/medical/med_007.htm"
    ]
    },
    {
    "cell_type": "code",
    "execution_count": 18,
    "metadata": {
    "collapsed": false
    },
    "outputs": [
    {
    "name": "stdout",
    "output_type": "stream",
    "text": [
    "2.30600413503\n"
    ]
    }
    ],
    "source": [
    "t = stats.t.isf(alpha / 2, n) # 両側検定を行う場合は, 2で割らないといけない点に注意\n",
    "print(t)"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "### スミルノフグラブス検定の有意点\n",
    "\n",
    "[スミルノフグラブス検定の有意点の表(片側)](http://aoki2.si.gunma-u.ac.jp/lecture/Grubbs/Grubbs-table.html)\n",
    "\n",
    "両側の場合は、危険率を2で割って表と照らし合わせる必要がある。\n",
    "データ数が10で、危険率5%でスミルノフグラブス検定を行う際の有意点は表より、 **2.127** になる。\n",
    "\n",
    "有意点を求める方法は、 [wikipedia](https://ja.wikipedia.org/wiki/%E5%A4%96%E3%82%8C%E5%80%A4) の通り。\n",
    "\n",
    "> 標本数を n、所要の有意水準を α、自由度 n - 2 のt分布の α / n × 100 パーセンタイルを t として、$\\frac{N-1}{\\sqrt{N}} \\times \\sqrt{\\frac{t^2}{N-2 + (t^2)}}$ を有意点とする。"
    ]
    },
    {
    "cell_type": "code",
    "execution_count": 22,
    "metadata": {
    "collapsed": false
    },
    "outputs": [
    {
    "name": "stdout",
    "output_type": "stream",
    "text": [
    "2.126645087195628\n"
    ]
    }
    ],
    "source": [
    "from scipy import stats\n",
    "from math import sqrt\n",
    "\n",
    "t = stats.t.isf((alpha/n) / 2, n-2) # 両側なので、 2で割る\n",
    "Gtest = (n-1)/sqrt(n) * sqrt(t**2 / (n-2 + t**2))\n",
    "print(Gtest)"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "OK! ここまで理解すれば後は簡単。\n",
    "\n",
    "データの平均値から最も離れた値に対して、平均を引いて標準偏差で割った値$\\frac{x - \\mu}{\\sigma}$を統計検定量として有意点と比較し、 `統計検定量 > 有意点` なら外れ値とみなす。これを再帰的に行い、外れ値がなくなるまで繰り返す。"
    ]
    },
    {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
    "## 参考リンク\n",
    "\n",
    "- https://ja.wikipedia.org/wiki/%E5%A4%96%E3%82%8C%E5%80%A4\n",
    "- https://sites.google.com/site/scriptofbioinformatics/mian-qiangmemo/waire-zhi-jian-chu-zhi-shi\n",
    "- http://www3.ocn.ne.jp/~stat/medical/med_007.htm"
    ]
    }
    ],
    "metadata": {
    "kernelspec": {
    "display_name": "Python 3",
    "language": "python",
    "name": "python3"
    },
    "language_info": {
    "codemirror_mode": {
    "name": "ipython",
    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.4.3"
    }
    },
    "nbformat": 4,
    "nbformat_minor": 0
    }