12_一般の分布に関する検定法.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 母比率の検定\n",
    "\n",
    "試行回数$n$、成功確率$\\theta$の二項分布に従う確率変数$X$を考える。  \n",
    "いま、与えられた比率$\\theta_0 \\in (0,1)$に対して、$\\theta=\\theta_0$を帰無仮説とする仮説検定を行いたいとする。  \n",
    "ここで$\\hat{\\theta}=X/n$は$\\theta$の最尤推定量である。  \n",
    "ここで棄却限界値をいかに定めるかが問題となるが、これには最尤推定量の漸近正規性を利用することが多い。  \n",
    "すなわち、帰無仮説のもとで$n \\rightarrow \\infty$のとき\n",
    "$$\n",
    "\\frac{\\sqrt{n} (\\hat{\\theta} - \\theta_0)}{\\sqrt{\\theta_0(1-\\theta_0)}}\n",
    "$$\n",
    "が標準正規分布に分布収束することを利用する。  \n",
    "ここで分母にある$\\theta_0(1-\\theta_0)$は$\\sqrt{n} (\\hat{\\theta} - \\theta_0)$の分散である。    \n",
    "すると、近似的に有意水準$\\alpha$の検定方式は、\n",
    "$$\n",
    "\\frac{\\sqrt{n} (\\hat{\\theta} - \\theta_0)}{\\sqrt{\\theta_0(1-\\theta_0)}} \\geq z_{\\alpha}\n",
    "$$\n",
    "とすれば得られる。ここで$a_{\\alpha}$は標準正規分布の上側$100 \\alpha %$点である。  \n",
    "両側対立仮説$\\theta \\neq \\theta_0$の場合、上と同様にして有意水準$\\alpha$の棄却域は\n",
    "$$\n",
    "\\left|\\frac{\\sqrt{n} (\\hat{\\theta} - \\theta_0)}{\\sqrt{\\theta_0(1-\\theta_0)}} \\right| \\geq z_{\\alpha /2}\n",
    "$$\n",
    "となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "統計量Z:2.219\n",
      "棄却限界値:1.960\n",
      "H0を棄却しH1を採択する\n"
     ]
    }
   ],
   "source": [
    "# サンプルデータで信頼度95%で両側検定を実行\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "theta0 = 0.5 # 母比率\n",
    "\n",
    "sample_n = 1000\n",
    "sample_x = 535\n",
    "theta_hat = sample_x/sample_n\n",
    "\n",
    "z_val = (theta_hat-theta0)/np.sqrt((theta_hat*(1-theta_hat))/sample_n)\n",
    "print(f\"統計量Z:{z_val:.3f}\")\n",
    "\n",
    "alpha = 0.05\n",
    "critical_val = abs(stats.norm.ppf(alpha/2))\n",
    "print(f\"棄却限界値:{critical_val:.3f}\")\n",
    "if(-critical_val<=z_val and critical_val>=z_val):\n",
    "    print(\"H0を棄却しない\")\n",
    "else:\n",
    "    print(\"H0を棄却しH1を採択する\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 母比率の差の検定\n",
    "\n",
    "試行回数$n_1$、成功確率$\\theta_1$の二項分布に従う確率変数$X_1$と、試行回数$n_2$、成功確率$\\theta_2$の二項分布に従う確率変数$X_2$が観測され、帰無仮説$\\theta_1=\\theta_2$の検定を行いたいとする。  \n",
    "つまり母比率に関する2標本問題である。ここで帰無仮説において$\\theta_1=(\\theta2)$の値そのものは決めていないことに注意する。  \n",
    "帰無仮説に限らず$\\hat{\\theta_1}-\\hat{\\theta_2}$の平均は$\\theta_1-\\theta_2$、分散は$\\theta_1(1-\\theta_1)/n_1 + \\theta_2(1-\\theta_2)/n_2$であるから、中心極限定理より$n_1$,$n_2$がともに大きければ\n",
    "$$\n",
    "\\frac{(\\hat{\\theta_1}-\\hat{\\theta_2})-(\\theta_1-\\theta_2)}{\\sqrt{\\frac{\\theta_1(1-\\theta_1)}{n_1} + \\frac{\\theta_2(1-\\theta_2)}{n_2}}} \n",
    "$$\n",
    "は近似的に標準正規分布に従う。  \n",
    "ここで分母に現れる$\\theta_1$,$\\theta_2$を$\\hat{\\theta_1}$,$\\hat{\\theta_2}$に置き換え、また帰無仮説のもとで$\\theta_1-\\theta_2$=0に代入することにより、有意水準$\\alpha$の検定方式\n",
    "$$\n",
    "\\frac{(\\hat{\\theta_1}-\\hat{\\theta_2})-(\\theta_1-\\theta_2)}{\\sqrt{\\frac{\\hat{\\theta_1}(1-\\hat{\\theta_1})}{n_1} + \\frac{\\hat{\\theta_2}(1-\\hat{\\theta_2})}{n_2}}} \\geq z_{\\alpha}\n",
    "$$\n",
    "が得られる。両側検定の場合も同様である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "統計量Z:3.816\n",
      "棄却限界値:1.960\n",
      "H0を棄却しH1を採択する\n"
     ]
    }
   ],
   "source": [
    "# サンプルデータで信頼度95%で両側検定を実行\n",
    "a_n = 1000\n",
    "a_x = 535\n",
    "theta_a = a_x/a_n\n",
    "\n",
    "b_n = 1000\n",
    "b_x = 450\n",
    "theta_b = b_x/b_n\n",
    "\n",
    "z_val = (theta_a-theta_b)/np.sqrt(theta_a*(1-theta_a)/a_n+theta_b*(1-theta_b)/b_n)\n",
    "print(f\"統計量Z:{z_val:.3f}\")\n",
    "\n",
    "alpha = 0.05\n",
    "critical_val = abs(stats.norm.ppf(alpha/2))\n",
    "print(f\"棄却限界値:{critical_val:.3f}\")\n",
    "if(-critical_val<=z_val and critical_val>=z_val):\n",
    "    print(\"H0を棄却しない\")\n",
    "else:\n",
    "    print(\"H0を棄却しH1を採択する\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ポアソン分布に関する検定\n",
    "\n",
    "パラメータ$\\lambda$のポアソン分布に従う確率変数$X$を考える。\n",
    "$$\n",
    "P(X=x) = \\frac{\\lambda^x}{x!} e^{- \\lambda}\n",
    "$$\n",
    "与えられた$\\lambda_0 > 0$に対して$\\lambda = \\lambda_0$を帰無仮説とする仮説検定を行いたいとする。  \n",
    "対立仮説が$\\lambda > \\lambda_0$の場合、棄却域を\n",
    "$$\n",
    "\\hat{\\lambda} - \\lambda_0 \\geq c\n",
    "$$\n",
    "と置くのが自然である。ここで$\\hat{\\lambda}=X$は$\\lambda$の最尤推定量である。  \n",
    "真の$\\lambda$の値が十分大きければ、$\\hat{\\lambda}$の分布は正規分布で近似される。  \n",
    "ポアソン分布の分散は$\\lambda$なので、帰無仮説の下では\n",
    "$$\n",
    "\\frac{\\hat{\\lambda} - \\lambda_0}{\\sqrt{\\lambda_0}}\n",
    "$$\n",
    "が近似的に標準積分布に従う。これに基づく棄却域は\n",
    "$$\n",
    "\\frac{\\hat{\\lambda} - \\lambda_0}{\\sqrt{\\lambda_0}} \\geq z_{\\alpha}\n",
    "$$\n",
    "である。両側検定も同様である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "統計量Z:2.582\n",
      "棄却限界値:1.960\n",
      "H0を棄却しH1を採択する\n"
     ]
    }
   ],
   "source": [
    "# サンプルデータで信頼度95%で両側検定を実行\n",
    "lambda0 = 15\n",
    "\n",
    "lambda_hat = 25\n",
    "\n",
    "z_val = (lambda_hat-lambda0)/np.sqrt(lambda0)\n",
    "print(f\"統計量Z:{z_val:.3f}\")\n",
    "\n",
    "alpha = 0.05\n",
    "critical_val = abs(stats.norm.ppf(alpha/2))\n",
    "print(f\"棄却限界値:{critical_val:.3f}\")\n",
    "if(-critical_val<=z_val and critical_val>=z_val):\n",
    "    print(\"H0を棄却しない\")\n",
    "else:\n",
    "    print(\"H0を棄却しH1を採択する\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 適合度検定\n",
    "\n",
    "カテゴリかるデータに関する検定には、尤度比検定のほかに**適合度検定**がよく用いられる。  \n",
    "カテゴリの個数$I$とし、それぞれの生起確率を$pi(1 \\leq i \\leq I)$とおく。  \n",
    "また総度数を$n$、カテゴリ$i$の観測度数を$x_i$とおき、帰無仮説のもとでの$p_i$の最尤推定量を$\\tilde{p_i}$とおく。  \n",
    "このとき、ピアソンのカイ二乗適合度検定量を\n",
    "$$\n",
    "T(X) = \\sum_{i=1}^{I} \\frac{(x_i - n \\tilde{p_i})^2}{n \\tilde{p_i}}, \\quad x=(x_1,...,x_I)\n",
    "$$\n",
    "と定義する。分母と分子に現れる$n \\tilde{p_i}$は期待度数とよばれる。  \n",
    "$T(X)$は帰無仮説の下で漸近的に自由度$d=I-1-k$のカイ二乗分布に収束する。  \n",
    "ここで$k$は帰無仮説のもとで自由に動けるパラメータの個数を表し、また$I-1$は生起確率に制約のないフルモデルにおいて自由に動けるパラメータの個数を意味する。  \n",
    "度数が十分に大きくない場合には、**イエーツの補正**が用いられることがある。  \n",
    "これは、$T(X)$の分子にある$(x_i - n \\tilde{p_i})^2$を\n",
    "$$\n",
    "(|x_i - n \\tilde{p_i}| - 0.5)^2\n",
    "$$\n",
    "に置き換える補正法である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "統計量T(x):9.458\n",
      "棄却限界値(右側):3.058\n",
      "棄却限界値(左側):0.465\n",
      "H0を棄却しH1を採択する\n",
      "chisquareの実行結果:Power_divergenceResult(statistic=9.458333333333334, pvalue=0.02377878015426507)\n"
     ]
    }
   ],
   "source": [
    "# サンプルデータで信頼度95%で両側検定を実行\n",
    "n = 100\n",
    "n_list = [55, 22, 16, 7] \n",
    "pi_hat_list = [0.4, 0.3, 0.2, 0.1] #A型:55、O型:22, B型:16, AB型:7\n",
    "\n",
    "chi_val = 0\n",
    "for i, item in enumerate(n_list):\n",
    "    chi_val += (item-100*pi_hat_list[i])**2/np.sum(100*pi_hat_list[i])\n",
    "print(f\"統計量T(x):{chi_val:.3f}\")\n",
    "\n",
    "alpha_l = 0.025\n",
    "alpha_r = 0.975\n",
    "critical_val_l = abs(stats.chi.ppf(alpha_l, len(n_list)-1))\n",
    "critical_val_r = abs(stats.chi.ppf(alpha_r, len(n_list)-1))\n",
    "print(f\"棄却限界値(右側):{critical_val_r:.3f}\")\n",
    "print(f\"棄却限界値(左側):{critical_val_l:.3f}\")\n",
    "if(chi_val<=critical_val_l or critical_val_r<=chi_val):\n",
    "    print(\"H0を棄却しH1を採択する\")\n",
    "else:\n",
    "    print(\"H0を棄却しない\")\n",
    "\n",
    "# 上記と同様の検定は、scipyのchisquareを使用すれば簡単に実行できる\n",
    "print(f\"chisquareの実行結果:{stats .chisquare(n_list, f_exp=[pi*n for pi in pi_hat_list])}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " 尤度比検定\n",
    "\n",
    "汎用的な検定の構成法として**尤度比検定**がある。  \n",
    "パラメータのベクトル$\\theta$が2つの部分に分かれており、$\\theta = (\\theta_1, \\theta_2)$とする。\n",
    "ただし、$\\theta_1$の次元(要素数)を$p$、$\\theta_2$の次元を$q$とする。  \n",
    "大きさ$n$の観測値$x_n$の確率密度関数を$f_n(x_n;\\theta_1, \\theta_2)$と書く。  \n",
    "帰無仮説を$H_0:\\theta_1=\\theta_{10}$とし、尤度比を\n",
    "$$\n",
    "\\lambda_n = \\frac{max_{\\theta_1,\\theta_2} f_n(x_n;\\theta_1, \\theta_2)}{max_{\\theta_2} f_n(x_n;\\theta_{10}, \\theta_2)}\n",
    "$$\n",
    "と定義する。  \n",
    "一般的に$H_0$のもとで$n \\rightarrow \\infty$のとき、$2 \\log \\lambda_n$の分布は自由度$p$のカイ二乗分布に従うことが示される。  \n",
    "これより$2 \\log \\lambda_n \\geq \\chi^2_{\\alpha}(p)$のときに$H_0$を棄却すれば、近似的に有意水準$\\alpha$となる。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.2 64-bit",
   "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.11.1"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "726ef1321fcd34e1bfe5e95996892c07db726849efbcd4abbb89e2457d1bbe6b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}