スミルノフグラブス検定

SmirnovGrubbsTest.ipynb

{
 "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
}

検定

• 推定: 母集団から抽出した標本をもとにして「母平均が◯◯と✕✕の間にある」というように未知の母数を予想
• 検定: 「母平均は△△である」という仮定をし、それが正しいか標本を元に判断する。

母集団の分布を仮定するわけだけど、これは「統計的仮説(Statistical Hypothesis)」という。
この仮説を採択(正しいとみなす)するか棄却(正しくないとみなす)するかを、標本を元に判断。
仮定が100%正しいわけではないので、危険率(=有意水準)という基準を設け、αで表す(通常は1%または5%)。

例えば、5%以下と少ない確率でしか起きないことが標本値をもとに計算した結果起きてしまった。
こうなったら仮説を棄却する。