Skip to content

Instantly share code, notes, and snippets.

@Tsutomu-KKE
Created July 4, 2014 22:29
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Tsutomu-KKE/a319bf6ceb52ac3700f0 to your computer and use it in GitHub Desktop.
Save Tsutomu-KKE/a319bf6ceb52ac3700f0 to your computer and use it in GitHub Desktop.
{
"metadata": {
"name": "",
"signature": "sha256:bb46f85200f566c8bc80f08923dafbcfa2241d3858b8035f59460ca07f521fb4"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#\u6b63\u898f\u5206\u5e03\u6027\u306e\u691c\u8a3c\n",
"\u6307\u5b9a\u3055\u308c\u305f\u5217\u306e\u6b63\u898f\u5206\u5e03\u3057\u3066\u3044\u308b\u5ea6\u5408\u3044\u3092$\\chi^2$\u691c\u5b9a\u3067\u8a08\u7b97\u3059\u308b\u3002 \u30ea\u30b9\u30c8\u306e\u500b\u6570\u306f\u300110\u306e\u500d\u6570\u3067\u300150\u4ee5\u4e0a\u3068\u3059\u308b\u3002\n",
"\n",
"\u6b63\u898f\u5206\u5e03\u306e\u9762\u7a4d\u304c10(=s)\u7b49\u5206\u306b\u306a\u308b\u3088\u3046\u306b\u5206\u5272\u3059\u308b\n",
"\n",
"\\begin{array}{rl}\n",
"n_s =& \\frac{n}{s} \\\\\n",
"Q =& \\sum_i{\\frac{(n_s - k_i)^2}{n_s}} \\\\\n",
"f(x, n) =& \\frac{1}{2^{\\frac{n}{2}}\\Gamma(\\frac{n}{2})}\n",
" \\int^x_0(u^{\\frac{n}{2} - 1} exp(-\\frac{u}{2}))du \\\\\n",
"p =& 1 - f(Q, s - 1)\n",
"\\end{array}\n",
"\n",
"N(0,1) \u306b\u304a\u3044\u3066\n",
"- P(x>0.253) = 0.4\n",
"- P(x>0.524) = 0.3\n",
"- P(x>0.842) = 0.2\n",
"- P(x>1.282) = 0.1"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#\u30c7\u30fc\u30bf\n",
"random.seed(1)\n",
"orig_data = random.randn(10000) * 10 + 50\n",
"orig_data = [sum(random.random(5)) for i in range(10000)]\n",
"hist(orig_data)\n",
"m, s = mean(orig_data), std(orig_data)\n",
"data = (orig_data - m) / s # \u6b63\u898f\u5316\n",
"cum = [len(data[data < i]) for i in (-1.282, -0.842, -0.524, -0.253, 0, 0.253, 0.524, 0.842, 1.282)]\n",
"num = [b - a for a, b in zip([0] + cum, cum + [len(data)])]\n",
"ns = len(data) / len(num)\n",
"Q = sum((i - ns)**2 / ns for i in num)\n",
"from scipy.integrate import quad\n",
"def f(x, n):\n",
" return quad(lambda u:u**(n / 2 - 1) * exp(-u / 2), 0, x)[0] / 2**(n / 2) / math.gamma(n / 2)\n",
"p = 1 - f(Q, len(num) - 1 - 2)\n",
"print('p = %.2f Q = %.2f %s' % (p, Q, num))\n",
"print('\u6b63\u898f\u5206\u5e03\u3067\u306a\u3044\u3068\u306f\u3044\u3048\u306a\u3044' if p > 0.05 else\n",
" '\u6b63\u898f\u5206\u5e03\u3067\u306a\u3044')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"p = 0.46 Q = 6.68 [1054, 972, 987, 991, 990, 984, 975, 991, 1035, 1021]\n",
"\u6b63\u898f\u5206\u5e03\u3067\u306a\u3044\u3068\u306f\u3044\u3048\u306a\u3044\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEMxJREFUeJzt3V2MnNV9x/HvjxhSQmwBSmsMuMKKjRJXSCAaiErTTN+Q\nU7VAb4BUSVGLolROQpqqVSEXZa/SpFJSyAXcQAKojVuUKAgSQjCEUYmq4hKZ4MS4vDRuWRebqNDY\nKILY4t+LeQyTZdld27szu3u+H2nFmTPnmXNmmPnN8ZnnJVWFJGl5O27cA5AkLTzDXpIaYNhLUgMM\ne0lqgGEvSQ0w7CWpATOGfZK1SR5K8sMkP0hyTVc/kWQyyfbu7wND21yX5Kkku5JcPFR/fpId3X03\nLtxTkiRNlZn2s09yGnBaVT2W5O3A94DLgMuBA1X1hSntNwJfAd4DnAE8AGyoqkqyDfh4VW1Lci/w\nxaq6b0GelSTp58w4s6+qvVX1WFd+CXiCQYgDZJpNLgW2VNXBqtoNPA1cmGQNsLKqtnXt7mDwpSFJ\nGoE5r9knOQs4D/i3ruoTSb6f5NYkJ3d1pwOTQ5tNMvhymFq/h9e/NCRJC2xOYd8t4XwV+GQ3w78Z\nWAecCzwHfH7BRihJOmYrZmuQ5Hjga8A/VNVdAFX1/ND9twD3dDf3AGuHNj+TwYx+T1cert8zta/1\n69fXM888c4RPQZKa90xVrZ+pwWx74wS4FdhZVTcM1a8ZavaHwI6ufDdwZZITkqwDNgDbqmovsD/J\nhd1jfhi46w2jfeYZqsq/Kq6//vqxj2Gx/Pla+Fr4Wsz8B7xztm+D2Wb2FwEfAh5Psr2r+zTwwSTn\nAgX8CPgoQFXtTHInsBM4BGyubiTAZuA24ETg3nJPHEkamRnDvqq+y/Sz/2/NsM1ngM9MU/894Jwj\nHaAk6dh5BO0i1ev1xj2ERcPX4nW+Fq/ztTgyMx5UNWpJajGNR5KWgiRU1XTHPr3Gmb0kNcCwl6QG\nGPaS1ADDXpIaYNhLUgNmPV2CtFitWnUqBw68OMIejwcOjqSnlStPYf/+F0bSl9rgrpdasgZn3hjl\n+2WU/QU/C5ord72UJAGGvSQ1wbCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCw\nl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJ\naoBhL0kNMOwlqQGGvSQ1wLCXpAbMGPZJ1iZ5KMkPk/wgyTVd/alJtiZ5Msn9SU4e2ua6JE8l2ZXk\n4qH685Ps6O67ceGekiRpqtlm9geBT1XVrwDvBT6W5N3AtcDWqjobeLC7TZKNwBXARmATcFOSdI91\nM3B1VW0ANiTZNO/PRpI0rRnDvqr2VtVjXfkl4AngDOAS4Pau2e3AZV35UmBLVR2sqt3A08CFSdYA\nK6tqW9fujqFtJEkLbM5r9knOAs4DHgFWV9W+7q59wOqufDowObTZJIMvh6n1e7p6SdIIrJhLoyRv\nB74GfLKqDry+MgNVVUlqvgY0MTHxWrnX69Hr9ebroSVpWej3+/T7/SPaJlUz53SS44FvAN+qqhu6\nul1Ar6r2dks0D1XVu5JcC1BVn+3a3QdcD/xX1+bdXf0HgfdX1Z9N6atmG4902GDSMcr3yyj7C34W\nNFdJqKrM1Ga2vXEC3ArsPBz0nbuBq7ryVcBdQ/VXJjkhyTpgA7CtqvYC+5Nc2D3mh4e2kSQtsBln\n9kl+HfgX4HFen9JcB2wD7gR+GdgNXF5V/9dt82ngT4FDDJZ9vt3Vnw/cBpwI3FtV10zTnzN7zZkz\ne2lgLjP7WZdxRsmw15Ew7KWBY17GkSQtD4a9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kN\nMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADD\nXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDVox7AJKms4IkI+tt\n5cpT2L//hZH1p9FLVY17DK9JUotpPDoyq1adyoEDL46411G+XzLC/kbZ16A/P3tLVxKqasbZgWGv\neTOYiS7X8B11f4a95m4uYe+avSQ1wLCXpAbMGvZJvpRkX5IdQ3UTSSaTbO/+PjB033VJnkqyK8nF\nQ/XnJ9nR3Xfj/D8VSdKbmcvM/svApil1BXyhqs7r/r4FkGQjcAWwsdvmpry+S8HNwNVVtQHYkGTq\nY0qSFsisYV9VDwPT7WIx3Y8BlwJbqupgVe0GngYuTLIGWFlV27p2dwCXHd2QJUlH6ljW7D+R5PtJ\nbk1ycld3OjA51GYSOGOa+j1dvSRpBI427G8G1gHnAs8Bn5+3EUmS5t1RHUFbVc8fLie5Bbinu7kH\nWDvU9EwGM/o9XXm4fs90jz0xMfFaudfr0ev1jmaIkrRs9ft9+v3+EW0zp4OqkpwF3FNV53S311TV\nc135U8B7quqPuh9ovwJcwGCZ5gFgfVVVkkeAa4BtwDeBL1bVfVP68aCqJcyDqpZqX4P+/OwtXXM5\nqGrWmX2SLcD7gXckeRa4HuglOZfBu/FHwEcBqmpnkjuBncAhYPNQem8GbgNOBO6dGvSSpIXj6RI0\nb5zZL9W+Bv352Vu6PF2CJAkw7CWpCYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwl\nqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIa\nYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwlqQGG\nvSQ1wLCXpAbMGvZJvpRkX5IdQ3WnJtma5Mkk9yc5eei+65I8lWRXkouH6s9PsqO778b5fyqSpDcz\nl5n9l4FNU+quBbZW1dnAg91tkmwErgA2dtvclCTdNjcDV1fVBmBDkqmPKUlaILOGfVU9DLw4pfoS\n4PaufDtwWVe+FNhSVQerajfwNHBhkjXAyqra1rW7Y2gbSdICO9o1+9VVta8r7wNWd+XTgcmhdpPA\nGdPU7+nqJUkjcMw/0FZVATUPY5EkLZAVR7ndviSnVdXebonm+a5+D7B2qN2ZDGb0e7rycP2e6R54\nYmLitXKv16PX6x3lECVpeer3+/T7/SPaJoOJ+SyNkrOAe6rqnO723wH/W1WfS3ItcHJVXdv9QPsV\n4AIGyzQPAOurqpI8AlwDbAO+CXyxqu6b0k/NZTxanAa/xY/y/99y7m/0z83P3tKVhKrKTG1mndkn\n2QK8H3hHkmeBvwE+C9yZ5GpgN3A5QFXtTHInsBM4BGweSu/NwG3AicC9U4NekrRw5jSzHxVn9kub\nM/ul2tegPz97S9dcZvYeQStJDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w\n7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDjvayhFoiVq06lQMHXhz3MCSNmRcvWeZGe0GR5Xwx\nkVH358VLNHdevESSBLiMIwmAFd2/Akdj5cpT2L//hZH1J8NeEgCHGOWy0YEDo/ti0YDLOJLUAMNe\nkhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWp\nAYa9JDXAsJekBhxT2CfZneTxJNuTbOvqTk2yNcmTSe5PcvJQ++uSPJVkV5KLj3XwkqS5OdaZfQG9\nqjqvqi7o6q4FtlbV2cCD3W2SbASuADYCm4CbkvgvC0kagfkI26nXF7sEuL0r3w5c1pUvBbZU1cGq\n2g08DVyAJGnBzcfM/oEkjyb5SFe3uqr2deV9wOqufDowObTtJHDGMfYvSZqDY73g+EVV9VySXwS2\nJtk1fGdVVZKZrmL8hvsmJiZeK/d6PXq93jEOUZKWl36/T7/fP6JtUjU/V5RPcj3wEvARBuv4e5Os\nAR6qqncluRagqj7btb8PuL6qHhl6jJqv8WggCdN8py5UbyPsa7n3t5yf26A/P+vzJwlVNXVJ/ecc\n9TJOkrclWdmVTwIuBnYAdwNXdc2uAu7qyncDVyY5Ick6YAOw7Wj7lyTN3bEs46wGvj6YObIC+Meq\nuj/Jo8CdSa4GdgOXA1TVziR3AjuBQ8Bmp/GSNBrztowzH1zGmX8u4yzV/pbzcxv052d9/izoMo4k\naekw7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLU\nAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNeBYLjiuo7Bq1akcOPDiuIchqTFecHzERnsBcPAi\n2Uu1v+X83Ab9LffP+ih5wXFJEuAyjqSxWNH9K3c0Vq48hf37XxhZf4uRYS9pDA4xymWjAwdG98Wy\nWLmMI0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNWCk58ZJ\nsgm4AXgLcEtVfW6U/U9VVTz88MO88sor4xyGJC24kZ3PPslbgP8AfgfYA/w78MGqemKozUjPZz85\nOcm6des56aT3jaS/V175b15++UnmdgKoPtCbh16XwznY+0z/Wiznc76/WV995ud9Mdf+Fsp89Ndn\n7q/F8QxOvrbwxnGGzbmcz36UM/sLgKerajdAkn8CLgWemGmjhfTqq6/y1rf+Ej/5ydYR9Xgj8Odz\nbNtnYT7US1EfX4vD+vhaHNZn7q/F6M6yuVjPsDnKNfszgGeHbk92dZKkBTbKmf2iuwbZcccdxyuv\n/JhVq/5gJP397Gf/ycsvj6QrSfo5o1yzfy8wUVWbutvXAa8O/0ib5GngnSMZkCQtH89U1fqZGowy\n7Fcw+IH2t4H/AbYx5QdaSdLCGNkyTlUdSvJx4NsMdr281aCXpNEY2cxekjQ+i+YI2iSbkuxK8lSS\nvx73eMYlyZeS7EuyY9xjGbcka5M8lOSHSX6Q5Jpxj2lckvxCkkeSPJZkZ5K/HfeYxinJW5JsT3LP\nuMcybkl2J3m8ez22vWm7xTCzn8sBV61I8j7gJeCOqjpn3OMZpySnAadV1WNJ3g58D7isxfcFQJK3\nVdVPu9+/vgv8ZVV9d9zjGockfwGcD6ysqkvGPZ5xSvIj4PyqmvFIrsUys3/tgKuqOggcPuCqOVX1\nMPDiuMexGFTV3qp6rCu/xOAAvNPHO6rxqaqfdsUTGPzuNdrDNBeJJGcCvwfcwuBQXM3hdVgsYe8B\nV5pRkrOA84BHxjuS8UlyXJLHgH3AQ1W1c9xjGpO/B/4KeHXcA1kkCnggyaNJPvJmjRZL2I9/LUmL\nVreE81Xgk90Mv0lV9WpVnQucCfxGkt6YhzRySX4feL6qtuOs/rCLquo84APAx7ql4DdYLGG/B1g7\ndHstg9m9GpfkeOBrwD9U1V3jHs9iUFU/Ab4J/Oq4xzIGvwZc0q1TbwF+K8kdYx7TWFXVc91/fwx8\nncGy+BsslrB/FNiQ5KwkJwBXAHePeUwasyQBbgV2VtUN4x7POCV5R5KTu/KJwO8C28c7qtGrqk9X\n1dqqWgdcCXynqv543OMalyRvS7KyK58EXAxMuyffogj7qjoEHD7gaifwzw3vcbEF+Ffg7CTPJvmT\ncY9pjC4CPgT8Zrdb2fbumggtWgN8p1uzfwS4p6oeHPOYFoPWl4BXAw8PvS++UVX3T9dwUex6KUla\nWItiZi9JWliGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDfh/CfLgHPSQBRYAAAAASUVO\nRK5CYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x6ae6588>"
]
}
],
"prompt_number": 2
}
],
"metadata": {}
}
]
}
@Tsutomu-KKE
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment