バートレット検定

BartlettTest.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from sklearn import datasets\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import bartlett"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "iris = datasets.load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "X = iris[\"data\"]\n",
    "ys = iris[\"target\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "data = [X[ys==0][:, 0], X[ys==1][:, 0], X[ys==2][:, 0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f2f4ec8a320>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE5NJREFUeJzt3X+wZ3V93/HnywUFMcAubIQAy5qEInGrFu6gRGJdiWlA\nC7ahFTJplG5mg41LTGbSMd1RiFObsU0nJTLDZgux0iaLlYqSFqg2rpFto9O7uIgGTRAhsEFdcF1E\n/LHgu398zx4vt/fHuXvvuefe3edj5jvf7znnc8/3zf2y39c95/M555OqQpIkgOcMXYAkaekwFCRJ\nLUNBktQyFCRJLUNBktQyFCRJLUNBktQyFCRJLUNBktQ6YugC5urEE0+stWvXDl2GJC0rO3fufKyq\nVs/WrtdQSPIbwK8ABdwLXFFV352w/XnATcA5wOPAm6rqwZn2uXbtWsbHx3urWZIORUke6tKut9NH\nSU4BrgLGqmodsAK4bFKzDcDeqvpJ4PeB9/ZVjyRpdn33KRwBHJ3kCOD5wN9O2n4J8IHm9S3ABUnS\nc02SpGn0FgpVtRv4PeBvgEeBfVX1sUnNTgEebto/DewDTuirJknSzPo8fbSS0ZHAi4AfA45J8ksH\nua+NScaTjO/Zs2chy5QkTdDn6aOfBb5SVXuqaj/wYeCnJ7XZDZwG0JxiOo5Rh/OzVNXWqhqrqrHV\nq2ftPJckHaQ+Q+FvgFcmeX7TT3ABcN+kNrcBb25eXwp8opz1R5IG02efwmcYdR7fzWg46nOArUne\nneTiptmNwAlJ7gd+E3hHX/VIkmaX5faH+djYWHmdgiTNTZKdVTU2W7tld0WzJB2MhRrtvtz+kJ4r\nQ0HSYWG2L/Mkh/wXfhfeEE+S1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIU\nJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1OotFJKcmWTXhMcTSd4+qc1rkuyb0OZd\nfdUjSZpdb9NxVtWXgJcDJFkB7AZunaLpXVX1hr7qkCR1t1injy4AvlxVDy3S+0mSDsJihcJlwLZp\ntp2X5J4kdyR5yVQNkmxMMp5kfM+ePf1VKUmHud5DIclzgYuBD02x+W7g9Kp6GfA+4CNT7aOqtlbV\nWFWNrV69ur9iJekwtxhHChcCd1fV1yZvqKonqurJ5vXtwJFJTlyEmiRJU1iMULicaU4dJTkpSZrX\n5zb1PL4INUmSptDb6COAJMcArwN+dcK6KwGqagtwKfDWJE8D3wEuq6rqsyZJ0vR6DYWq+jZwwqR1\nWya8vg64rs8aJEndeUWzJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiS\nWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWoaCJKllKEiSWr2FQpIzk+ya8HgiydsntUmSP0hy\nf5LPJTm7r3okSbPrbY7mqvoS8HKAJCuA3cCtk5pdCJzRPF4BXN88S5IGsFinjy4AvlxVD01afwlw\nU418Gjg+ycmLVJMkaZLFCoXLgG1TrD8FeHjC8iPNOknSAHoPhSTPBS4GPjSPfWxMMp5kfM+ePQtX\nnDQHSRbkIS1li3GkcCFwd1V9bYptu4HTJiyf2qx7lqraWlVjVTW2evXqnsqUZlZVsz66tJOWssUI\nhcuZ+tQRwG3ALzejkF4J7KuqRxehJknSFHobfQSQ5BjgdcCvTlh3JUBVbQFuBy4C7geeAq7osx5J\n0sx6DYWq+jZwwqR1Wya8LuDX+qxBktSdVzRLWvZWrVq1IAMA5ruPVatWDfybmL9ejxQkaTHs3bt3\nSXTiHwqjyzxSkCS1DAVJUstQkCS1DAVJUstQkCS1DAVJUstQkCS1DAVJUstQkCS1DAVJUstQkCS1\nDAVJUstQkCS1Zr1LapJXAdcApzftw2gqhB/vtzRJ0mLrcuvsG4HfAHYCz/RbjiRpSF1CYV9V3dF7\nJZKkwU3bp5Dk7CRnA9uT/Lsk5x1Y16yfVZLjk9yS5ItJ7kty3qTtr0myL8mu5vGuef73SJLmYaYj\nhX8/aXlswusCXtth/9cCd1bVpUmeCzx/ijZ3VdUbOuxLktSzaUOhqtYDJPnxqnpg4rYks3YyJzkO\neDXwlmZ/3we+P59iJUn96jIk9ZYp1n2ow8+9CNgDvD/JZ5PckOSYKdqdl+SeJHckeUmH/UqSejLt\nkUKSFwMvAY5L8o8nbDoWOKrjvs8GNlXVZ5JcC7wDeOeENncDp1fVk0kuAj4CnDFFLRuBjQBr1qzp\n8NaSDid19bFwzXFDlzGqY5lLVU29IbkEeCNwMXDbhE3fAm6uqv8z446Tk4BPV9XaZvlngHdU1etn\n+JkHgbGqemy6NmNjYzU+Pj7TW0uDScJ0/6bUn6Xye18qdUwlyc6qGput3Ux9Ch8FPprkvKr6i7kW\nUFVfTfJwkjOr6kvABcBfTiryJOBrVVVJzmV0Ouvxub6XJGlhdLlO4ReTXD5p3T5gvAmOmWwC/rgZ\nefQAcEWSKwGqagtwKfDWJE8D3wEuq6Uas5J0GOgSCs8DXswPO5d/AfgK8LIk66vq7dP9YFXt4tlD\nWQG2TNh+HXDdnCqWJPWmSyi8FHhVVT0DkOR64C7gfODeHmuTJC2yLkNSVwIvmLB8DLCqCYnv9VKV\nJGkQXY4U/i2wK8knGd0h9dXAv2muOfhfPdYmSVpks4ZCVd2Y5Hbg3GbVv6qqv21e/1ZvlUmSFl3X\nSXaew+jq5L3ATyZ5dX8lSZKG0mWSnfcCbwK+APygWV3Ap3qsS5I0gC59Cm8EzqwqO5V1SFu1ahV7\n9+6d936SHPTPrly5km984xvzrkE6WF1C4QHgSBxppEPc3r17B79FwXwCRVoIXULhKUajj/6MCcFQ\nVVf1VpUkaRBdQuE2nn1DPEnSIarLkNQPJDkaWNPc2E6SlpylcOpt5cqVQ5cwb7MOSU3yD4FdwJ3N\n8suTeOQgacmoqnk/FmI/h8IggS7XKVzD6MK1b0J7k7tZp+OUJC0/XUJhf1Xtm7TuB1O2lCQta106\nmr+Q5BeBFUnOAK4CZpx1TZK0PHU5UtjEaK7m7wHbgCeAaedQkCQtX11GHz0FbG4ekqRD2LShkORP\nGd3jaEpVdXEvFUmSBjPTkcLvLVoVkqQlYdpQqKo/n+/OkxwP3ACsY3TU8c+r6i8mbA9wLXARo9tp\nvKWq7p7v+0qSDk6X0UfzcS1wZ1VdmuS5wPMnbb8QOKN5vAK4vnmWJA2g6yQ7c5bkOEZTd94IUFXf\nr6pvTmp2CXBTjXwaOD7JyX3VJEmaWW+hALyI0Wxt70/y2SQ3NPM6T3QK8PCE5UeadZKkAfQ5+ugI\n4GxgU1V9Jsm1wDuAd861yCQbgY0Aa9asmeuPS5I66nP00SPAI1X1mWb5FkahMNFu4LQJy6c2656l\nqrYCWwHGxsaGnQVFkg5hvY0+qqqvJnk4yZnNLbcvAP5yUrPbgLcluZlRB/O+qnp0Pu8rSTp4s44+\nau539LvATwFHHVhfVV3ulLoJ+ONm5NEDwBVJrmx+fgtwO6PhqPczGpJ6xVz/AyRJC6fLkNT3A1cD\nvw+sZ/TF3amDurnN9tik1VsmbC/g1zpVKknqXZdQOLqq/ixJquoh4JokO4F39VybtKjq6mPhmuOG\nr0EaUJdQ+F6S5wB/neRtjDqCX9BvWdLiy+880c7ANVgNCXXNoCXoMNflNNCvM7oS+SrgHOCfAW/u\nsyhJ0jC63Dr7/wI0RwtXVdW3eq9KkjSIWY8UkowluRf4HHBvknuSnNN/aZKkxdalT+GPgH9RVXcB\nJDmf0Yikl/ZZmCRp8XUJhWcOBAJAVe1I8nSPNR3SRncLn5+hO0MlHbq6hMKfJ/lDRvMzF/Am4JNJ\nzgZw/oO5me0LPYlf+pIG0yUUXtY8Xz1p/d9jFBKvXdCKJEmD6TL6aP1iFCJJGl6X0UcvTHJjkjua\n5Z9KsqH/0iRJi63LxWv/CfifwI81y38FvL2vgiRJw+kSCidW1X8FfgBQVU8Dz/RalSRpEF1C4dtJ\nTqCZhS3JK4F9vVYlSRpEl9FHv8loMpyfSPK/gdXApb1WJUkaRJfRR3cn+fvAmUCAL1XV/t4rkyQt\nui6jj/4JozkVvgC8EfjggQvXJEmHli59Cu+sqm819zy6ALgRuL7fsqRhJBn0sXLlyqF/BYes2X73\nXdosxG1qlrpO9z5qnl8P/Meq+h9J/nWXnSd5EPhWs4+nq2ps0vbXAB8FvtKs+nBVvbvLvqWFthC3\nF/E2JUuXn0s3XUJhd3Pvo9cB703yPDrO0dxYX1WPzbD9rqp6wxz2J0nqSZcv93/K6OK1f1BV3wRW\nAb/Va1WSpEHMGgpV9VRVfbiq/rpZfrSqPtZx/wV8LMnOJBunaXNeM3HPHUle0nG/kqQedDl9NB/n\nV9XuJD8KfDzJF6vqUxO23w2cXlVPJrkI+AhwxuSdNIGyEWDNmjU9lyxJh6+59A3MWVXtbp6/DtwK\nnDtp+xNV9WTz+nbgyCQnTrGfrVU1VlVjq1ev7rNkSTqs9RYKSY5J8iMHXgM/B3x+UpuT0ozxSnJu\nU8/jfdUkSZpZn6ePXgjc2nznHwH8SVXdmeRKgKrawuh2GW9tpvf8DnBZOW5MkgbTWyhU1QP8cNa2\nieu3THh9HXBdXzVIkuam1z4FSdLyYihIklqGgiSpZShIklqGgiSpZShIklqGgiSpZShIklqGgiSp\nZShIklqGgiSpZSgsoFWrVs174naY/+Txq1atGvg3IWm56nuSncPK3r17l8Tk4AfCRZLmyiMFSVLL\nUJAktQwFSVLLUJAktQwFSVLLUJAktXoNhSQPJrk3ya4k41NsT5I/SHJ/ks8lObvPeiRJM1uM6xTW\nV9Vj02y7EDijebwCuL55liQNYOjTR5cAN9XIp4Hjk5w8cE2SdNjq+0ihgI8lKeAPq2rrpO2nAA9P\nWH6kWffoxEZJNgIbAdasWdNftfNUVx8L1xw3dBmjOiTpIPQdCudX1e4kPwp8PMkXq+pTc91JEyZb\nAcbGxoa/j8Q08jtPLJnbXNQ1Q1chaTnq9fRRVe1unr8O3AqcO6nJbuC0CcunNuskSQPoLRSSHJPk\nRw68Bn4O+PykZrcBv9yMQnolsK+qHkWSFsm2bdtYt24dK1asYN26dWzbtm3okgbV5+mjFwK3Nnfs\nPAL4k6q6M8mVAFW1BbgduAi4H3gKuKLHeiTpWbZt28bmzZu58cYbOf/889mxYwcbNmwA4PLLLx+4\numFkKZwDn4uxsbEaH///LnlYEpIsnT6FJVDH4cjf/fKybt063ve+97F+/fp23fbt29m0aROf//zk\nExvLW5KdVTU2a7vl9j+wobB86jjULNQ8FX42S8eKFSv47ne/y5FHHtmu279/P0cddRTPPPPMgJUt\nvK6hMPR1CtKyUVUL8tDScdZZZ7Fjx45nrduxYwdnnXXWQBUNz1CQdNjavHkzGzZsYPv27ezfv5/t\n27ezYcMGNm/ePHRpg3E6TkmHrQOdyZs2beK+++7jrLPO4j3vec9h28kM9iksqKVyLn+p1CFp6eja\np+CRwgJbqM7I+Vi5cuXQJUhapgyFBbQQf537V76kIdnRLElqGQqSpJahIElqGQqSpJahIElqGQqS\npJahIElqGQqSpJahIElqGQqSpJahIElq9R4KSVYk+WyS/z7Ftrck2ZNkV/P4lb7rkSRNbzFuiPfr\nwH3AsdNs/2BVvW0R6pAkzaLXI4UkpwKvB27o830kSQuj79NH/wH4l8APZmjzC0k+l+SWJKf1XI8k\naQa9hUKSNwBfr6qdMzT7U2BtVb0U+DjwgWn2tTHJeJLxPXv29FCtJAn6PVJ4FXBxkgeBm4HXJvkv\nExtU1eNV9b1m8QbgnKl2VFVbq2qsqsZWr17dY8mSdHjrLRSq6rer6tSqWgtcBnyiqn5pYpskJ09Y\nvJhRh7QkaSCLPh1nkncD41V1G3BVkouBp4FvAG9Z7HokST+U5TYf8NjYWI2Pjw9dRm+co1lSH5Ls\nrKqx2dp5RbMkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJai36dwuEuybzbOGRVUl8MhUXmF7qkpczT\nR5KklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWotu0l2kuwBHhq6jh6dCDw2dBE6aH5+\ny9eh/tmdXlWzTnK/7ELhUJdkvMvsSFqa/PyWLz+7EU8fSZJahoIkqWUoLD1bhy5A8+Lnt3z52WGf\ngiRpAo8UJEktQ2EJSfLGJJXkxUPXou6SPJNkV5J7ktyd5KeHrkndJTkpyc1JvpxkZ5Lbk/ydoesa\niqGwtFwO7GietXx8p6peXlUvA34b+N2hC1I3GU1zeCvwyar6iao6h9Fn+MJhKxuOobBEJHkBcD6w\nAbhs4HJ08I4F9g5dhDpbD+yvqi0HVlTVPVV114A1DcrpOJeOS4A7q+qvkjye5Jyq2jl0Uerk6CS7\ngKOAk4HXDlyPulsH+O9sAo8Ulo7LgZub1zfjKaTl5MDpoxcDPw/c1JyWkJYdh6QuAUlWAY8Ae4AC\nVjTPp5cf0JKX5MmqesGE5a8Bf7eqvj5gWeogyQXA1VX16qFrWSo8UlgaLgX+c1WdXlVrq+o04CvA\nzwxcl+aoGTm2Anh86FrUySeA5yXZeGBFkpcmOWz/7RkKS8PljEZATPTf8BTScnF0MyR1F/BB4M1V\n9czQRWl2zZH4PwJ+thmS+gVGo8e+Omxlw/H0kSSp5ZGCJKllKEiSWoaCJKllKEiSWoaCJKllKEiS\nWoaCJKllKEiSWv8PnxtPP+qg4OoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2f51109a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot(data)\n",
    "plt.xticks([1, 2, 3], [\"A\", \"B\", \"C\"])\n",
    "plt.ylabel(\"sepal length\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "k = len(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "N_i = [d.shape[0] for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sum_N_i_minus_1 = sum([N_i[i] - 1 for i in range(k)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "C = 1 + 1 / (3 * (k - 1)) * (sum([1 / (N_i[i] - 1) for i in range(k)]) - 1 / sum_N_i_minus_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "s_i = [np.std(d, ddof=1) for d in data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "T = sum([(N_i[i] - 1) * np.log(sum([(N_i[i] - 1) * (s_i[i] ** 2) for i in range(k)]) / sum_N_i_minus_1) for i in range(k)])\n",
    "T -= sum([(N_i[i] - 1) * np.log(s_i[i] ** 2) for i in range(k)])\n",
    "T /= C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.005701874401502"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BartlettResult(statistic=16.005701874401502, pvalue=0.00033450760701630349)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bartlett(data[0], data[1], data[2])"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}