Skip to content

Instantly share code, notes, and snippets.

@Wybxc

Wybxc/方差分析.ipynb

Created June 1, 2023 15:18
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save Wybxc/d7973b7b5442ef09c2efe5b45843cb74 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save Wybxc/d7973b7b5442ef09c2efe5b45843cb74 to your computer and use it in GitHub Desktop.
Download ZIP
双因子方差分析计算器
Display the source blob
Display the rendered blob
Raw
方差分析.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[[15., 15., 17.],\n",
" [19., 19., 16.],\n",
" [16., 18., 21.]],\n",
"\n",
" [[17., 17., 17.],\n",
" [15., 15., 15.],\n",
" [19., 22., 22.]],\n",
"\n",
" [[15., 17., 16.],\n",
" [18., 17., 16.],\n",
" [18., 18., 18.]],\n",
"\n",
" [[18., 20., 22.],\n",
" [15., 16., 17.],\n",
" [17., 17., 17.]]])"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"\n",
"data_raw = \"\"\"\n",
"15，15，17 19，19，16 16，18，21\n",
"17，17，17 15，15，15 19，22，22\n",
"15，17，16 18，17，16 18，18，18\n",
"18，20，22 15，16，17 17，17，17\n",
"\"\"\".strip()\n",
"\n",
"data = np.array([[s.split(\"，\") for s in l.split()] for l in data_raw.splitlines()])\n",
"data = np.float64(data)\n",
"data\n"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(4, 3, 3)"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S, T, R = data.shape\n",
"S, T, R"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([17.33333333, 17.66666667, 17. , 17.66666667]),\n",
" array([17.16666667, 16.5 , 18.58333333]),\n",
" array([[15.66666667, 18. , 18.33333333],\n",
" [17. , 15. , 21. ],\n",
" [16. , 17. , 18. ],\n",
" [20. , 16. , 17. ]]),\n",
" 17.416666666666668)"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mu_a = data.mean(axis=(1, 2))\n",
"mu_b = data.mean(axis=(0, 2))\n",
"mu_r = data.mean(axis=2)\n",
"mu = data.mean(axis=(0, 1, 2))\n",
"\n",
"mu_a, mu_b, mu_r, mu\n"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sa=2.7500000000000124, sb=27.16666666666663, sr=41.33333333333333, sab=73.50000000000004, s=144.75\n"
]
}
],
"source": [
"sa = ((mu_a - mu) ** 2).sum() * T * R\n",
"sb = ((mu_b - mu) ** 2).sum() * S * R\n",
"sr = ((data - mu_r[:, :, np.newaxis]) ** 2).sum()\n",
"s = ((data - mu) ** 2).sum()\n",
"sab = s - sa - sb - sr\n",
"print(f\"sa={sa}, sb={sb}, sr={sr}, sab={sab}, s={s}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"na=3, nb=2, nr=24, nab=6, n=35\n"
]
}
],
"source": [
"na = S - 1\n",
"nb = T - 1\n",
"nab = (S - 1) * (T - 1)\n",
"nr = S * T * R - 1 - na - nb - nab\n",
"print(f\"na={na}, nb={nb}, nr={nr}, nab={nab}, n={na+nb+nr+nab}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fa=0.5322580645161314, fb=7.887096774193538, fab=7.112903225806456\n"
]
}
],
"source": [
"if R == 1:\n",
" fa = (sa / na) / (sab / nab)\n",
" fb = (sb / nb) / (sab / nab)\n",
" fab = 0\n",
" print(f\"fa={fa}, fb={fb}\")\n",
"else:\n",
" fa = (sa / na) / (sr / nr)\n",
" fb = (sb / nb) / (sr / nr)\n",
" fab = (sab / nab) / (sr / nr)\n",
" print(f\"fa={fa}, fb={fb}, fab={fab}\")"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pa=0.6645281856425167, pb=0.002329798210769174, pab=0.00019217323641351758\n"
]
}
],
"source": [
"from scipy import stats\n",
"\n",
"if R == 1:\n",
" pa = 1 - stats.f.cdf(fa, na, nab)\n",
" pb = 1 - stats.f.cdf(fb, nb, nab)\n",
" print(f\"pa={pa}, pb={pb}\")\n",
"else:\n",
" pa = 1 - stats.f.cdf(fa, na, nr)\n",
" pb = 1 - stats.f.cdf(fb, nb, nr)\n",
" pab = 1 - stats.f.cdf(fab, nab, nr)\n",
" print(f\"pa={pa}, pb={pb}, pab={pab}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<thead>\n",
"<tr><th>方差来源 </th><th style=\"text-align: right;\"> 平方和</th><th style=\"text-align: right;\"> 自由度</th><th>F值 </th><th>显著性 </th></tr>\n",
"</thead>\n",
"<tbody>\n",
"<tr><td>因子A </td><td style=\"text-align: right;\"> 2.75 </td><td style=\"text-align: right;\"> 3</td><td>F(3, 24)=0.532</td><td>p=0.665 </td></tr>\n",
"<tr><td>因子B </td><td style=\"text-align: right;\"> 27.1667</td><td style=\"text-align: right;\"> 2</td><td>F(2, 24)=7.887</td><td>p=0.002 </td></tr>\n",
"<tr><td>交互作用 </td><td style=\"text-align: right;\"> 73.5 </td><td style=\"text-align: right;\"> 6</td><td>F(6, 24)=7.113</td><td>p=0.000 </td></tr>\n",
"<tr><td>误差 </td><td style=\"text-align: right;\"> 41.3333</td><td style=\"text-align: right;\"> 24</td><td> </td><td> </td></tr>\n",
"<tr><td>总和 </td><td style=\"text-align: right;\">144.75 </td><td style=\"text-align: right;\"> 35</td><td> </td><td> </td></tr>\n",
"</tbody>\n",
"</table>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from IPython.display import HTML, display\n",
"import tabulate\n",
"\n",
"headers = [\"方差来源\", \"平方和\", \"自由度\", \"F值\", \"显著性\"]\n",
"if R == 1:\n",
" table = [\n",
" [\"因子A\", sa, na, f\"F({na}, {nab})={fa:.3f}\", f\"p={pa:.3f}\"],\n",
" [\"因子B\", sb, nb, f\"F({nb}, {nab})={fb:.3f}\", f\"p={pb:.3f}\"],\n",
" [\"误差\", sab, nab, \"\", \"\"],\n",
" [\"总和\", s, na + nb + nab, \"\", \"\"],\n",
" ]\n",
"else:\n",
" table = [\n",
" [\"因子A\", sa, na, f\"F({na}, {nr})={fa:.3f}\", f\"p={pa:.3f}\"],\n",
" [\"因子B\", sb, nb, f\"F({nb}, {nr})={fb:.3f}\", f\"p={pb:.3f}\"],\n",
" [\"交互作用\", sab, nab, f\"F({nab}, {nr})={fab:.3f}\", f\"p={pab:.3f}\"],\n",
" [\"误差\", sr, nr, \"\", \"\"],\n",
" [\"总和\", s, na + nb + nab + nr, \"\", \"\"],\n",
" ]\n",
"display(HTML(tabulate.tabulate(table, headers=headers, tablefmt=\"html\")))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "base",
"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.10.9"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment