Chapter4ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 4 Analysis of Variance\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 One-way ANOVA\n",
    "\n",
    "### 4.1.1 A motivating example: the Spock trial\n",
    "\n",
    "In 1968 Dr. Benjamin Spock was tried in Boston for conspiring against the government for helping young men to escape the military draft. He was convicted by the Boston federal court, but the judgement was overturned by the Court of Appeals in 1969 for many reasons, one of which was cited as the bias of the presiding judge Francis Ford. \n",
    "\n",
    "Dr. Spock, a pediatrician, was very famous for his books on rearing of children, and thus was widely admired by women. As a matter of fact, the jury in Spock trial has no women. Note that jury panels, though randomly selected, should reflect the demographics. \n",
    "\n",
    "In any particular trial, there may not be any woman on the jury, but it is worthwhile to examine if the jury panels of Judge Ford had fewer women than other judges in Boston in few months before the trial. Data are available for jury panels for 7, but we investigate the data for only 4 judges including Judge Ford. \n",
    "\n",
    "We can start our analysis with a visualization of the data set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDT09PZ2dnh4eHp6enw8PD////iz9LxAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAacklEQVR4nO3d4WKa2rqGUWI0Js1OYu7/ZnfVxpiuNhV4gQ8Y48de\n5Bx1TpGnCqI270BvzdQTgCUQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFB\ngJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAk\nCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKE\nBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGA\nkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQI\nEBIECAkChAQBQoIAIUGAkCBASBAgJAgYIaQGZqbDVp4PZ4IhIElIECAkCBASBAgJAoQEAUKC\nACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgpMo6fYKZKQippLu792NGH/9DeUKq\n6O7uWNJ5Lax9XcyEkAq6O/lYCStfGTMhpIKEND9Cquj00k5IcyKkks4HG06Lq18X8yCkuhy1\nmxEhVeZ9pNkQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQI\nEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQE\nAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQ\nIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBo4b08rhrjnb7l6GGgEmM\nGNLhvvm0HWQImMiIIe2bzY/X09Lb86bZDzEETGTEkDbN62X5tdkMMQRMZMSQmuZvf8SGgIl4\nRoKAcfeRnt9OS/aRWJoxD39vr47a3R8GGQKmMe77SPvT+0ib3aP3kVgWZzZAQJ2QmmvDDAFD\nmSSkf4YiJGZGSBAw6huyN796ExIzM2JILxshsVRjvrQ77Jrt6R1ZL+1YmnH3kX40zY93IbE8\nIx9seNs2u4OQWJzRj9o9NptnIbE04x/+fr3/9xuuQmJmpngf6UFILE2dU4RGHgKShAQBQoIA\nIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\nICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJ\nAoTEkt3djTSQkFiwu7uxShISy3V3N1pJQmLGmu+dQ/ruErmZjHKVgkOwBl7aQUDjYAP0N96G\nJCQWTEgwK0KCACFBgJDmYrTjT0tiH4nfjPeOyJIIia9GPNllSYS0EP84heV2/z7ZZeBTYOZJ\nSAvR/C/l2FH7ay1mRXYkpIXIhfS/Lh2tPqTxCGlQwZA6WcyKLE9IgxLSWghpUEKaln2khRDS\ntIS0EEKalpAWQkjTEtJCCGlaQloIIa2FkAYlpLUQ0qCEtBZCGpSQpmUfaSGENC0hLYSQpiWk\nhRDStIS0EEKalpAWQkhrIaRBCWkthDQoIa2FkAa17JDu7qp/2559pIVYdEh39b8jTEgLseSQ\n7u7qlySkhZh7SP/+qr2xflqy4/wrjySk2809pG/Vf0IS0vBDjGPRIf0sabyflqxOSINadkgL\neqB6E9KghLQWQhrU0kOqzj7SQghpWkJaCCFNS0gLsfSQqj9QQloIIU1LSAshpLUQ0qCEtBZC\nGpSQ1kJIg1p6SNXZR1qIjiF1+r1YIf2XkBaiW0jdfsFcSP8lpIXoFNL50wlVQ/pyvnf1B0pI\nC7G8kL5+Aqn6AyWkhVjcS7vfPsu3mAeqNyENqvrBhu8+Jf7N58tbX63qR9FzhDSo6oe/288v\n92x5y/zmQ0iDyv3TPcy/+B1CT3a0oDeMhTSoqTsaIKQoIQ1LSELKrP+Bb7/XSEK6XZd9kDE3\nVCFNOJKQbtd6Q83uy48bUoeZD/ZA/zpEL6SFaLuhBt+LvWlDTYbUZeZDPdDjf3GlkAa1opA6\nTX2gB3qCr4AV0qBW9NJOSCNcpeAQ41jTwYaKL+3sIy1E9aNiDjakCGlQqwppgPn1Xv8D336v\nkYR0u+obavX59V7/A99+r5GEdLvqG2r1+fVe/wPffq+RhHS76htq9fnNh5AGVX1DrT6/+RDS\noKpvqNXnNx9CGlT1DbX6/Hqv/4Fvv9dIQrpd9Q21+vx6r/+Bb7/XSEK6XfUNtfr8eq//gW+/\n10hCul31DbX6/Hqv/4Fvv9dIQrpd9Q21+vx6r/+Bb7/XSEK6XfUNtfr85kNIg6q+oVaf33wI\naVDVN9Tq85sPIQ2q+oZafX691//At99rpO6Te3ncnb4jard/GWqIYqpvqP3n1++DiELq4HB/\n9X1r20GGKGfxIfX8aLyQOtg3mx+vp6W3502zH2KIcsqH1FPv79Qfev0PfPu9Ruo6uU3zell+\nbTZDDFGOkISUvMr5es3f/ogNUc7SQ2r6/sjLKI/CGDwjDap8SL1HqH2wYTzj7iM9v52W7CON\nZe2Hv8cz5uHv7dVT+v1hkCGqqb6hVp9f7/U/8O33GqnH+0j70/tIm92j95FqbKjV59d7/Q98\n+71GcmbD7apvqNXn13v9D3z7vUYaaHKLPJhTfUOtPr/e6//LX0N+HfgEIT1tmvunYYcoo/qG\nWn1+vdf/9R+DfrH+mCG97prN0/vj6SnHKUIlNtTq80sa9icqRgzp9VTQvnk4vL/tmm+fk4Q0\n0oZafn5Bnc7CuP2RvvmSfa5y8nB872h/fif20NwPMUQ55TfU5c/v8x3jQX+/afRThJrd1R/p\nIcpZwYZafH7X9Qz5+019Q3q8v/lJ8HyRH+fXdE4RGsXqQ+r5W6KjhfTY4tXkw3Hv6Ozw4BSh\nUQhpJiFtvj9o8MVhc6mt+f4JSUgpqw+p5ycPRwup3Vun+498Nt8+HwkpRkj9Tk8fLaRd8+3J\np10JKURIw87vaqY3X/KPV3nbbP9x/mknQhppQzC/fvO7munNl/zjVbq8ddVyiFmrviGYX7/5\nXc305kv+8SpC+l7gqNOgG0L1DbX6/K5mevMl+1yl4BDj6Lsh9Py2q9lvqNXndzXTmy/Z5yoF\nhxhHzw2h57sg899Qq8/vaqY3X/IvV3neHV/V7d7a387NQ8xYp7MrP/X+tqu5b6jV53c105sv\n+eerbM+PVrOJlrSYkPoa9CM07/U31Orzu5pph5V/tfzUbA/HkJ6ah/Y3dNsQq9YM21H5DbX6\n/K5m2mHlXy1vmsP57AZH7QYx+CdIi2+o1ed3NdMOK/96uXkX0oCEVHt+VzPtsPKvlu9/PSO9\nfv9BvT5DMKDqG2r1+V3NtMPKv1r+tY/03OYs8JZDMKDqG2r1+V3NtMPKv/5j9+sw6/dfZtJr\nCIZTfUOtPr+rmXZY+V/+Or6P1Ox+tL+Z24dYMftIted3NdMOK7/9VQoOMQ9Cqj2/q5l2WPnt\nr1JwiHkQUu35Xc20w8pvf5WCQ8yDkGrP72qmHVb+9R/7zY1ndXUfYs2EVHt+VzPtsPKvlvc3\nnx7ZeQgGVH1DrT6/q5l2WPnXy9n3j/40BAOqvqFWn9/VTDus/OvlYTZ5IY2k+oZafX5XM+2w\n8q+W975FaFD2kWrP72qmHVb+9R/bbfYjfX8YYs2EVHt+VzPtsPKv/3h2sGFI9UI6fvK95zey\ntNhQVxNSm+/+7jjEqpUL6e6u9/dItNlQVxNS+KzvPw2xatVCurvLliSkj2VH7Wat7XepfIbU\n9pp/8a/5rSWkR0ft1uWjo5GGW01I74+++3tdjg2N1tF6Qrr9SbrzEKtWfUVU24cT0t+HWLWr\nFTHe00ALQrrMtMOda3+VgkPMw+eKGHHHpAUhXWba4c61v0rBIebhsiJG3cW/nZAuM+1w5778\n9eP4pcW+s6GjlsedWx8+Hnz+Q9/+akLa/npAfYvQwEo+IQ2uX0j93zYeLaSnZvP88z++1254\na+yoX0iBEzBGC+m+eT391zetMoQO50r889VwK7fPtMOd+3JHf1+IENJMDL6PtJaQPp+RNu1v\n6LYhKKxySE3ilMDbZ9rhzl0t20dat3JH7b7sFc3oYIOjdutWLaTgR6XGDenl/Yfv/l6xYiFF\nP3Q4akjN5tF3NjCY1YT0cHxV9yP+mSQhcbKal3a/ThF6eG5/M7cPwWr1O9gwq5De398e75tm\ns29/Q7cPQVHF9pHiRj77+/DQ5oh7pyEoSUiXmXa4c1//fD0+ITXbx/Y3dPMQVCWky0w73Lmr\n5efjz7rc78O7SEKaCyFdZtrhzl0vN83utf1NtBmCwoR0mWmHO3e1/HzcO/r5jBQ+AC4kTlYT\n0k8vpx/tu3fUjrw1hfTTi6N2DGJVIR2Oh+3uHbVbI/tIl5l2uHNf/jqd2bAPf92qkGZCSJeZ\ndrhzV8unc+3SB7+FNBtCusy0w527Xnb296oJ6TLTDnfuanmIb9B/F9JsCOky0w53rv1VCg7B\nHKwvJIe/GYCQ+hESJ+sLKUtIM2Ef6TLTDneu/VUKDkGCkC4z7XDn2l+l4BAkCOky0w537vqP\nw/74Faubffb0byHNhJAuM+1w566W3zanowxNs4m+LyukmRDSZaYd7tzV8rZ5OD4XHfbNrv0N\n3TYEK7aakPwaBUNaTUib5rxzdBASA1hNSPtmezzd7mXb+ITsGtlHusy0w527/sOvUayakC4z\n7XDnvvx1+jWKbfTXkYQ0G0K6zLTDnWt/lYJDkCCky0w73Lmr5V101+iPQ1CYkC4z7XDnrpeH\nWZNC4mQ1Id038d9G+n0IVmw1IR122yE+bS4kTlYTUpcfUm85BIXZR7rMtMOdu14W0qoJ6TLT\nDneu/VUKDkGCkC4z7XDn2l+l4BAkCOky0w537stfz7vjq7pd9msihTQTcwipz+8zj3yu3c//\nmw/2MYBASHd3PUoaLaSnZnv6BMVT89D+hm4bgnm4uxvgRvuHdHfXp6RRP4/068Pm7W/otiGY\nhePmmr/VprdzSJ2vfvtMO9y5L3f0XUgrdnmgztvrlFP5m5HmFThF6NjQa3Mfm9K7kGajekjH\n+Y0zrcw+0vOmiX4iafiQCj7oc/T5QJXsaMR/kfsetdv9ei05s0/I1nzY5+fqgSq5QucT0ul9\npGb3IzSdPw6RV/SFyPxUfw0+o5AGERui96Gc1ESYhJAG9vGEJBQyMl9+8hD+PebRDjYIiYzU\n13FFv7F4vO1bSP1Yfx96f0Hk5vhkNL/D32MPtFDV199s9pE2zevpv96QXafqD9RsQvIl+utW\n/YGaTUj7yzOSn3VZo+oP1GxCen887SO9bGZ2ZsPoAzGJ2YQ00NuXQmJmhAQBKz2zYfSBFsr6\n+yAkeqi+/mazjzSQ6o8Pv1R/oITELFR/oITELFR/oIQ0kuobAv0IaSRCIkNIECAkerD+PgiJ\nHqqvP/tIzEL1B0pIzEL1B0pIzEL1B0pII6m+IdCPkEYiJDKEBAFCogfr74OQ6OFz/dX8TQL7\nSMzC5YEq+useQmIWPh6oqr+TIyRKuOkHi7/7nZyp5195pCWFNPUDvRjFnpA+JiOkkQgppVhH\no09HSCzOFHtso4b08nj+7ebd/mWoIVoS0hItPKTD/dVO6PffFS4k+rh0tMh9pH2z+XH+7Yq3\n502zH2IIOFn0wYaPHyU7em02QwwBXywypC9vKnz/DoOQiFhkSBWfkRS7bIsM6fjDzW+npTr7\nSEIiY8zD39uro3b3h0GGaEtIZIz7PtL+9D7SZvfofSSWxZkNLNgi95H+dbOFThpmIZYZ0uGh\nabbPv27E4W+Gt8iQDpvziXbnGxESw1tkSPvm6WdNT5vTaXZFQlLssi0ypM35im+b+zchsTAT\nnCJ02G6FxMKMGNJ98/Em7P1WSCzLiCE9NQ+/lt6arZAYwSL3kd73l3qe//FWke2biGWG9P66\n+1h6exASw1toSJWGYA2ENBLFLpuQRiIkMoQEAUKCACGxYPaRIEBIECCkEZX6HQWihDSSptwv\n+5AkpJE0ZX+0kXkRkpAIWH1IXtqRICQHGxbMPhIECAkChAQBQhqJYpdNSCMREhlCggAhQYCQ\nTryXtEz2kcbl7IaFEtKonG+3VEK69YIR55A6XXXItUBvQrr1gv+LOHbU6YpCqk1It14wE9L/\nOnYkJH4RUi9C4kxIQiJASEJaMPtIt15QSHxDSLdeUEh8Q0i3XlBIfENIt15QSHxDSLdeUEiU\nICQhESAkIREgJCEtmH2kWy8oJH5z/YEYId16QSHx1ZePlgnp1gsKiS++fkhTSLdeUEir1vtD\nmrmZjHKVwYYQEr+Z6FsDhCSkhZnm2zeEJCQChCQkAoQkJAKEJCQChCQkAoQkJAKEJCQChCQk\nAoQkJAKEJCQChCQkAoQkJAKEJCQChNT5J12ExCchdf6RMSHxafUhnT9RKST6WXFI53yERMJ6\nQ/rox0s7AuYeUmefX5LR8ffMw1+ewbwJqZ8h1wIzstqQmkhHQuJs7iF130fq9f6RfSS+WnFI\nCULiTEhCIkBIQiJASEIiQEhCIkBIQiJASEIiQEhCIkBIQiJASEIiQEhCIkBIQiJASEIiQEhC\nIkBIQiJgtSElPo0kJD6sNaReX3kiJH630pD6fQmXkPidkIREwEpD8tKOrLWGFCpJSJytNaTQ\nazshcTb3kPp8q53vtSNm5iF1d+xolIFYhdWG9H7sSEiErDekcQdi4YQEASsPCTKEBAFCgoCV\nh6RYMoQEAUKCACFBgJAgYOUhQYaQIEBIELDykBRLhpAgQEgQICQIEBIErDwkyBASBAgJAlYe\nkmLJEBIECAkChAQBQoKAlYcEGUKCACFBwMpDUiwZQoIAIUGAkCBASBCw8pAgQ0gQMGpIL4+7\n5mi3fxlqCJjEiCEd7ptP20GGaE2xZIwY0r7Z/Hg9Lb09b5r9EEO0JiQyRgxp07xell+bzRBD\ntCYkMkYMqWn+9kdsiNaERIZnJAgYdx/p+e20VGcfCTLGPPy9vTpqd38YZAiYxrjvI+1P7yNt\ndo/eR2JZVn5mg2LJqBNSc22YIf4w6FgDsXB1Qhp5iJEHYuGEBAFCgoBRz2y4eTfI9s3MjBjS\nk5BYrDFf2r1uvv/wRGAImMao+0iv358YlBiiJcWSMe7Bhqer81YHGqIdIZHhqB0ECAkChAQB\nKw8JMoQEAUKCgJWHpFgyVhvS3d1IA7EKaw3p7u5UkpDIWGlId3fnkoRExsJDav7iHNLf/r8T\nfOadmVt4SH/166UdZKw1pHcdkbTakCBJSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBCw7pB8Sz4hKw7p7u44\njJRIWG9IH79HISQCVhvSx0+NKYkEIQmJgNWG9PlTY0Kiv/WG9K4jclYc0vmAnaN2JKw5JO8j\nEbPukCBESBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\nUDQkmJkOW3k+nMlUvy/m10/p+ZWeXEvV74v59VN6fqUn11L1+2J+/ZSeX+nJtVT9vphfP6Xn\nV3pyLVW/L+bXT+n5lZ5cS9Xvi/n1U3p+pSfXUvX7Yn79lJ5f6cm1VP2+mF8/pedXenItVb8v\n5tdP6fmVnlxL1e+L+fVTen6lJ9dS9ftifv2Unl/pybVU/b6YXz+l51d6cjAXQoIAIUGAkCBA\nSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAhYU0lPp+/J032z2h6ln8VeHh6Z5\neJ16Ft97qfwAV55bO69dfkNgNPvTjxxsypa0Oc2vdEmHTeUHuPLcWnndVA7ptXk4HJ8zH6ae\nyF/sjzPbN7up5/GdXeUHeDEhPTXbyut5d55b2SlumuNzZdnpHf3o9LNFo6k8tzaafe3N4Kz4\nFJvN1DP4u7fa/1IuJqTX8lvpT4dmO/UUvrNvnqaewt9tm7fSD3DlubVUej2fPDXPU0/h736+\ndNpPPYe/e2x+1H6AK8+tpdLr+ehtU3ln/mm3aR6nnsTfvB6Pg5R+gCvPraXS6/n9ePi29Au7\nnx7Kvra7P75xUPoBrjy3lkqv55+291PP4F8OVY82PJxeE5d+gCvPraXS6/n97X77NvUc/qnq\nKmwupp7JX9WdWWuF1/L7+3PtA3bn95HemqJPmkIaU+G1fHwXZOopfOt0ZsNhV3Yf6aTyAyyk\ncTxU/xf1fK5d7drrrr13IY2k/kuT/aa5L/18VPsBXlJIMB0hQYCQIEBIECAkCBASBAgJAoQE\nAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQ\nIEBIECAkCBASBAgJAoQEAUKCACHNwp9/rK70T9itjIdiFoRUnYdiFoRUnYdiFoRUnYdiFk7J\nnLs5/+9+0+z/u/h032yK/zT5UglpFn4Padv8tPtcfDgt7o6LzXbKia6WkGbht5B+NJvX99fN\ncfH5enF7eD9sm+eJJ7tKQpqF30LaNS/vx3BOi8+fi4efi4dmN+VM10pIs/BbSL+OMvxn8ZfJ\nprliVvosCKk6K7263aWhf4c01RwRUnnNcc/nM6SXzx2jl//sIznMMBkhVbdr9u+PpwMI983T\n8aDc74fqLounY3nvTw42TEFI1b2cdnuOR+mePt882n2+ebT93DE6L27epp3wOgmpvOf75v7l\ntPS4aR5+7Qk9fjmzYftyObOhedDRFIS0DM5nmJiQZq5pfry/H447UkxJSDP3eN5D2kw9j7UT\n0tw9bZvm3vPR1IQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBwP8B5YfShqbjOQ0AAAAASUVORK5C\nYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Spock <- read.csv(file=\"../Data/SpockTrial.csv\", header=TRUE, sep=\",\")\n",
    "Spock$Judge<-as.factor(Spock$Judge);\n",
    "# Box plot with jittered points \n",
    "boxplot(perc.women~Judge,data=Spock)\n",
    "stripchart(perc.women~Judge, vertical = TRUE, data = Spock,\n",
    "    method = \"jitter\", add = TRUE, pch = 20, col = 'blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.2 One-way ANOVA\n",
    "\n",
    "In the Spock trial data, let $r=4$ denote the number of judges,  $Y_{ij}$ be the percentage of women in the $j$th panel for the $i$th judge. Let Judge 1 be the judge in the Spock trial. We can propose the following model, for $j=1,\\ldots, n_i, i=1,\\ldots, r$,\n",
    "$$\n",
    "Y_{ij}=\\mu_i + \\epsilon_{ij}, \n",
    "$$\n",
    "where $\\{\\epsilon_{ij}\\}$ are i.i.d. $N(0,\\sigma^2)$. In the Spock trial data, $r=4$, $n_2=6$, and $n_1=n_3=n_4=9$. This model is a one-way **an**alysis **o**f **va**riance model in its cell means form.  We will discuss other forms later in this chapter.  \n",
    "\n",
    "[\n",
    "\n",
    "🔴IMPORTANT❗🔴\n",
    "\n",
    "$Y_{ij}$: (observed) outcome\n",
    "\n",
    "$\\epsilon_{ij}$: (unobserved) error term\n",
    "\n",
    "$\\mu_i$: cell mean -> pupulation mean of cell i]\n",
    "\n",
    "e.g.:\n",
    "\n",
    "cell 1 $i=1: Y_{1,1},\\cdots,Y_{1,5}$ we have cell mean $\\mu_1$ and $n_1=5$\n",
    "\n",
    "cell 2 $i=2: Y_{2,1},\\cdots,Y_{2,19}$ we have cell mean $\\mu_2$ and $n_2=19$\n",
    "\n",
    "cell 3 $i=3: Y_{3,1},\\cdots,Y_{3,10000}$ we have cell mean $\\mu_3$ and $n_3=10000$\n",
    "\n",
    "cell 4 $i=4: Y_{4,1},\\cdots,Y_{4,12}$ we have cell mean $\\mu_4$ and $n_4=12$\n",
    "\n",
    "]\n",
    "\n",
    "Because $n_1,n_2, n_3$ and $n_4$ are not equal, this is an imbalanced ANOVA model.🔴IMPORTANT❗🔴 If $n_1=n_2=\\cdots=n_r$, then the ANOVA model is **balanced**.\n",
    "\n",
    "In this model, $\\mu_i$ represents the mean percentage of women in the panels for the $i$th judge, and $\\sigma^2$ represents the variance in the percentages across panels. It is easy to see that, by using one $\\sigma^2$ across all judges, we assume the panels of all four judges have the same amount of variability. \n",
    "\n",
    "The question of interest in the Spock trial data can now be translated to whether $\\{\\mu_i\\}_{i=1}^r$ are the same, where $\\{\\mu_i\\}$ and $\\sigma^2$ are unknown. \n",
    "\n",
    "\n",
    "The estimators for $\\mu_i, i=1,2,\\ldots, r$ are simply the within-group sample means, i.e., for $i=1,\\ldots, r$,\n",
    "$$\n",
    "\\hat{\\mu}_i = \\bar{Y}_{i\\cdot} = \\frac{1}{n_i}\\sum_{j=1}^{n_i} Y_{ij}.\n",
    "$$\n",
    "We have two observations on our estimators $\\hat{\\mu}_i$, $i=1,\\ldots, r$. \n",
    "1. The estimator $\\hat{\\mu}_i$ is also the maximum likelihood estimator for $\\mu_i$. \n",
    "2. $\\hat{\\mu}_i$ is the best linear unbiased estimator if $\\{Y_{ij}\\}$ are mutually uncorrelated but not necessarily normally distributed.\n",
    "\n",
    "We can estimate the variance of errors as \n",
    "$$\\hat{\\sigma}^2 \\equiv \\frac{1}{n_T-r} \\sum_{i=1}^{r} \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\hat{\\mu}_i\\big)^2, $$\n",
    "\n",
    "[\n",
    "\n",
    "🔴IMPORTANT❗🔴\n",
    "\n",
    "$Y_{ij}-\\hat{\\mu}_i$ is residual\n",
    "\n",
    "]\n",
    "\n",
    "\n",
    "where $n_T\\equiv \\sum_{i=1}^r n_i$ is the (total) sample size.\n",
    "\n",
    "We can call the `aov()` function to fit a one-way ANOVA model in `R`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value   Pr(>F)    \n",
       "Judge        3 1591.3   530.4   17.61 1.06e-06 ***\n",
       "Residuals   29  873.5    30.1                     \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "anova.fit<- aov(perc.women~Judge,data=Spock)\n",
    "\n",
    "# Summary\n",
    "summary(anova.fit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.3 Alternative forms of ANOVA\n",
    "\n",
    "**Factor effect form**\n",
    "\n",
    "$\\sum_{i=1}^r w_i=1. w_{i}>0\\ \\forall i$.\n",
    "\n",
    "Let $\\mu=\\sum_{i=1}^r w_i \\mu_i$ and $\\tau_i = \\mu_i -\\mu$. Then $\\sum_{i=1}^r w_i \\tau_i=0$. We can rewrite the ANOVA model as \n",
    "$$\n",
    "Y_{ij}=\\mu+\\tau_i+\\epsilon_{ij}, \n",
    "$$\n",
    "where $\\{\\epsilon_{ij}\\}$ are i.i.d. $N(0,\\sigma^2)$. \n",
    "\n",
    "The least squares estimators are $\\hat{\\mu}=\\sum_{i=1}^r w_i \\bar{Y}_{i\\cdot}$ and $\\hat{\\tau}_i = \\bar{Y}_{i\\cdot} -\\hat{\\mu}$. \n",
    "\n",
    "[\n",
    "$w_i$ 就是根据全局分布对所得数据进行靠近全局分布的改良\n",
    "]\n",
    "\n",
    "Notes in `R`:\n",
    "1. One of $\\tau_i$ is set to zero by default.\n",
    "2. `R` assumes $w_i=n_i/n_T$  for one-way ANOVA, but equal weights for higher-order ANOVAs. \n",
    "    \n",
    "**Regression form**\n",
    "  \n",
    "There are multiple equivalent forms to turn a cell-mean model into the typical linear regression form. \n",
    "In the Spock trial data, the regression equation takes the following form, for $j=1,\\ldots, n_i,\\ \\quad i=1,\\ldots, 4,$ \n",
    "$$\n",
    "  Y_{ij}=\\mu +\\tau_1 X_{1,ij}+\\tau_2 X_{2,ij} + \\tau_3 X_{3,ij}  +\\epsilon_{ij}. \n",
    "$$\n",
    "where $\\{\\epsilon_{ij}\\}$ are i.i.d. $N(0,\\sigma^2)$. Here we set $\\tau_4$ to be zero, and thus we ignore the the fourth term $ X_{4,ij}$.\n",
    "\n",
    "\n",
    "There are multiple choices in the coding of $\\{X_{1,ij},\\ldots, X_{4,ij}\\}$. \n",
    "1. *Dummy variables.* $X_{l,ij}=1$ when $l=i$, and 0 otherwise. \n",
    "2. *Unweighted effect coding.* For $l=1,2,3$,  $X_{l,ij}=1$ when $i=l$, $X_{l,ij}=-1$ when $i=4$, and $X_{l,ij}=0$ otherwise. This is equivalent to an ANOVA model with equal weights. \n",
    "3. *Weighted effect coding.* For $l=1,2,3$, $X_{l,ij}=1$ when $i=l$, $X_{l,ij}=-n_l/n_4$ when $i=4$, and $X_{l,ij}=0$ otherwise. This is equivalent to an ANOVA model with unequal weights. \n",
    "\n",
    "It is easy to see that all three regression models are equivalent, since they are all equivalent forms of the same ANOVA model. However, the interpretations of the model parameters (i.e., $\\tau$s) may differ slightly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = perc.women ~ as.factor(Judge), data = Spock)\n",
       "\n",
       "Residuals:\n",
       "    Min      1Q  Median      3Q     Max \n",
       "-10.300  -3.300  -0.100   3.078  11.983 \n",
       "\n",
       "Coefficients:\n",
       "                  Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept)         14.622      1.829   7.993 8.16e-09 ***\n",
       "as.factor(Judge)2   18.994      2.893   6.567 3.41e-07 ***\n",
       "as.factor(Judge)3   14.478      2.587   5.596 4.85e-06 ***\n",
       "as.factor(Judge)4   12.178      2.587   4.707 5.73e-05 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 5.488 on 29 degrees of freedom\n",
       "Multiple R-squared:  0.6456,\tAdjusted R-squared:  0.6089 \n",
       "F-statistic: 17.61 on 3 and 29 DF,  p-value: 1.057e-06\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lm.fit<- lm(perc.women~as.factor(Judge),data=Spock)\n",
    "summary(lm.fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.4 Sum of squares decomposition \n",
    "\n",
    "Define the residuals $e_{ij}=Y_{ij}-\\hat{Y}_{ij}=Y_{ij}-\\bar{Y}_{i\\cdot}$, $j=1,2,\\ldots, n_i$ and $i=1,2, \\ldots, r$. \n",
    "From the model equation, we can see that $e_{ij}$ is an estimator of $\\epsilon_{ij}$.\n",
    "\n",
    "A simple identity holds that \n",
    "$$\n",
    "Y_{ij}-\\bar{Y}_{\\cdot\\cdot} = (\\bar{Y}_{i\\cdot}-\\bar{Y}_{\\cdot \\cdot})+(Y_{ij}-\\bar{Y}_{i\\cdot}).\n",
    "$$\n",
    "Squaring both sides of the equation yields \n",
    "\\begin{align}\n",
    "\\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\bar{Y}_{\\cdot \\cdot}\\big)^2 = & \\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(\\bar{Y}_{i\\cdot}-\\bar{Y}_{\\cdot \\cdot}\\big)^2 +\\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\bar{Y}_{i \\cdot}\\big)^2\\\\\n",
    "=  & \\sum_{i=1}^r n_i \\big(\\bar{Y}_{i\\cdot}-\\bar{Y}_{\\cdot \\cdot}\\big)^2 +\\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\bar{Y}_{i \\cdot}\\big)^2,\n",
    "\\end{align}\n",
    "where the first equality holds because $2\\sum_{i=1}^r\\sum_{j=1}^{n_i} \\big(\\bar{Y}_{i\\cdot}-\\bar{Y}_{\\cdot \\cdot}\\big)\\big(Y_{ij}-\\bar{Y}_{i \\cdot}\\big)=0$. \n",
    "\n",
    "Some terminologies:\n",
    "\n",
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "    <th> </th>\n",
    "    <th>Acronym</th>\n",
    "    <th>Definition</th>\n",
    "     <th>d.f.</th>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Total sum of squares </td>\n",
    "    <td>SSTO</td>\n",
    "    <td>$\\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\bar{Y}_{\\cdot \\cdot}\\big)^2$</td>\n",
    "    <td> $n_T-1$ </td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td> Treatment sum of squares | Across group variation</td>\n",
    "    <td>SSTR </td>\n",
    "    <td >$\\sum_{i=1}^r n_i \\big(\\bar{Y}_{i\\cdot}-\\bar{Y}_{\\cdot \\cdot}\\big)^2$  </td>\n",
    "    <td>  $r-1$ </td>\n",
    "  </tr>\n",
    "    <tr>\n",
    "    <td> Residual sum of squares |Within group variation </td>\n",
    "    <td>  SSE  </td>\n",
    "    <td> $\\sum_{i=1}^r \\sum_{j=1}^{n_i} \\big(Y_{ij}-\\bar{Y}_{i \\cdot}\\big)^2$ </td>\n",
    "    <td>   $n_T-r$ </td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n",
    "We can see that \n",
    "$$\n",
    "{\\rm SSTO} \\ = \\ {\\rm SSTR} \\ +\\ {\\rm SSE}, \\quad {\\rm df(SSTO)}\\ = \\ {\\rm df(SSTR)} \\ + \\ {\\rm df(SSE)}.\n",
    "$$\n",
    "\n",
    "We can define the mean squares accordingly, \n",
    "\\begin{align}\n",
    "{\\rm MSTR}\\ = & \\ {\\rm SSTR}/{\\rm df(SSTR)} \\ = &  \\ {\\rm SSTR}/(r-1)\\\\\n",
    "{\\rm MSE}\\ = & \\ {\\rm SSE}/{\\rm df(SSE)} \\ = &  \\ {\\rm SSE}/(n_T-r)\\\\\n",
    "{\\rm MSTO}\\ = & \\ {\\rm SSTO}/{\\rm df(SSTO)} \\ = &  \\ {\\rm SSTO}/(n_T-1)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value   Pr(>F)    \n",
       "Judge        3 1591.3   530.4   17.61 1.06e-06 ***\n",
       "Residuals   29  873.5    30.1                     \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Summary\n",
    "summary(anova.fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.5 Gauss-Markov properties\n",
    "\n",
    "We assume that $Y=X\\beta + \\epsilon$, where $\\epsilon$ is an $n \\times 1$ vector of mean zero, uncorrelated errors with a common variance $\\sigma^2$. Assume that $X \\in \\mathbb{R}^{n \\times r}$ is of rank $r$ and $n>r$. The least squares estimator of $\\beta$ is $\\hat{\\beta} = (X^T X)^{-1}X^T Y$. Consequently, $\\hat{Y}=X\\hat{\\beta}$ and $e=Y-X\\hat{\\beta}$. \n",
    "\n",
    "Properties:\n",
    "1. $\\mathbb{E}\\big[\\hat{Y} \\big]=X{\\beta}$ and ${\\rm cov}\\big(\\hat{Y}\\big)= X(X^T X)^{-1} X^T \\sigma^2  \\equiv H \\sigma^2$.\n",
    "2. $\\mathbb{E}[e]=0$ and ${\\rm cov}(e)=(I-H)\\sigma^2$.\n",
    "3. ${\\rm cov}\\big(\\hat{Y},e\\big)=0$. Furthermore,  if $\\epsilon$ is i.i.d. $N(0,\\sigma^2)$, $\\hat{Y}$ and $e$ are independent. \n",
    "4. $\\mathbb{E}[{\\rm SSE}]=(n_T-r)\\sigma^2$, where ${\\rm SSE}=e^T e$. Thus, $\\mathbb{E}[{\\rm MSE}]=\\sigma^2$.\n",
    "\n",
    "How does the ANOVA model relate to the Gauss-Markov models? \n",
    "- $Y$: turn $\\{Y_{ij}\\}$ to $Y \\in \\mathbb{R}^{n_T}$\n",
    "\n",
    "- $X$: an $n_T \\times r$ matrix of $0$'s and $1$'s\n",
    "\n",
    "- $\\beta$: $\\beta\\equiv (\\mu_,\\mu_2,\\ldots, \\mu_r)^T$\n",
    "\n",
    "Corresponding properties of ANOVA:\n",
    "1. $\\big\\{ \\hat{Y}_{ij}\\big\\}$ are independent of the residuals $\\{e_{ij}\\}$\n",
    "2. ${\\rm MSTR} \\perp {\\rm MSE}$\n",
    "3. $\\mathbb{E}[{\\rm MSE}]=\\sigma^2$\n",
    "4. $\\mathbb{E}[{\\rm MSTR}]=\\sigma^2 + \\sum n_i (\\mu_i-\\mu)^2/(r-1)$ where $\\mu=\\sum_{i=1}^r (n_i/n_T)\\mu_i$\n",
    "5. ${\\rm SSE}/\\sigma^2 \\sim \\chi^2_{n_T-r}$ and when $\\mu_1 = \\cdots = \\mu_r$ and ${\\rm SSTR}/\\sigma^2 \\sim \\chi^2_{r-1}$.\n",
    "6. When $\\mu_1 = \\cdots = \\mu_r$ then $F={\\rm MSTR}/{\\rm MSE}$ has an $F$-distribution with degrees of freedom $(r-1,n_T-r$).  \n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.6 F-test\n",
    "\n",
    "Recall that the question of interest in the Spock trial data can now be translated to whether $\\{\\mu_i\\}_{i=1}^r$ are the same. Formally, we want to test the null hypothesis \n",
    "$$ H_0: \\mu_1 = \\mu_2 = \\cdots = \\mu_r $$\n",
    "against the alternative \n",
    "$$ H_1: {\\rm not \\ all\\ } \\mu_i\\ {\\rm are\\ the\\ same}.$$\n",
    "\n",
    "【\n",
    "\n",
    "Full model:\n",
    "$I:Y_{ij}=\\mu_i+\\epsilon_{ij},i=1,\\cdots,r,j=1,\\cdots,r_i$\n",
    "\n",
    "Reduced model:\n",
    "$II: Y_{ij}+\\mu+\\epsilon_{ij}$\n",
    "\n",
    "】\n",
    "\n",
    "The first test statistic we consider is the F-statistic $F^*\\equiv  {\\rm MSTR}/{\\rm MSE}$. Under the null, $F^*$ follows an F-distribution with d.f.s $(r-1,n_T-r)$ when assuming $\\epsilon$s are i.i.d. $N(0,\\sigma^2)$. We reject the null hypothesis at the significance level $\\alpha$, if $F^* > F(1-\\alpha; r-1, n_T-r)$. Here the value of $F(1-\\alpha; r-1, n_T-r)$ is known as the *critical value*, and the set $(F(1-\\alpha; r-1, n_T-r), +\\infty)$ is the *rejection region*.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD/pQD///+afi2oAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAgAElEQVR4nO3dh5aiSBiG4SIY2sCO93+zC5gwE/6vqKLe5+yZddqW\nom3eIYjoTgAmc3PPALAEhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJ\nMEBIgAFCAgwQEmCAkAADHkJysPXff//NPQtLN2Iptw9nhiHSUoc09ywsHCElgZDUCCkJhKRG\nSEkgJDVCSgIhqRFSEghJjZCSQEhqhJQEQlIjpCQQkprPkI4rl21Op23usrVoCLxHSGoeQ6qy\n5kyK7aY9oaKQDIEPCEnNY0hrV6+H1plbVaeqvW0/BD4gJDWPIWXtA52r2v9liiHwASGpeQzJ\nufufp+9nyxKSMUJSm2GN1PxZsUbyipDUZthHWleX2/ZD4ANCUuOoXRIISY3XkZJASGqc2ZAE\nQlIjpCQQkhohJYGQ1OYKideRvCIktXBCmnhtI3yjDcnrZa+CxKZdIsQhCacdBUJKBSFJEVIq\nCEnKd0j7TdluUZbrvWoIvEVIUn5DqvLO3hmnCHlFSFJ+Q1q77O/Q3jruMk5a9YqQpPyGlLnD\n7faBt1F4RUhSfkN6eHWIF2S9IiQp1kipICQp7/tIu2N7i30k3whJyvPh76Jz1C6vJEPgPUKS\n8v460rp9HSkrN7yO5BchSXFmQyoISYqQUkFIUoSUCkKSIqRUEJIUIaWCkKQIKRWEJEVIqYgr\npMMqc6vdwwjNELvbrbcz8f2ks1/z2F5tsfvi5jpzxXkeqpVzq8P7h12mTkiJiCqk9eU1+2Nn\nhHqI3F1vvZ+JSSGdx8zuJZ3PHtg0N9srBLtvJRFSKmIKaeOyek1Q1f87Pnz9RwuTQjo0n9p1\n2rrV9QtbV1TNqujQNLZq/ii/TZ2QEhFRSMdrQKv7Yn0eRRhS6Z6+q3D7dl7WzQqp+jUBQkpF\nRCGtz9tT9Tqp3DbLb5XXK4N6Mb5e9qodrNmBaXPble5yJfnrkl65vP1/Xi//j/d2P6Cr2SXa\nvv4g7unW/Z3cX9+uQEip8B7Sf/28mVjxsDfiXN3C+jmk4rpDc/5ok/atBLcICtckdqwbeLq3\nE1L59noHVaebp1XU2r1m9/AEEFISfIfUs6N3c/W4DeWafZXnDv6aL66aQJz7a/7quo/7a1dp\nG7d7vvc+gV27B1S4hyODzY7R7Qt5W+P+MtE/9/NtP4SUhIjWSM8h7a9fvHdQNl/sfOjjY0in\ndtsudy/3difQ1Fk9HUA4Zve/b1xZnQ7FZaLbMrtucH56AggpCRHtIz2HdPvfvYPutxx3m+Ip\npFW9NjmeVyEP93YncLvU8F2VdTf12mPeZXeiX7btCCkVEYVU3vaRdtXvkIpbEPcv7uu1x7pd\nkz3e+yOkIu/+rWo/Fe/+HV8/9ZiQUhFRSJvrRtS+2UT7EdLK5dvd8SmkU5Y3/73c+36VdnHM\ni+PLFw/uHte349+ElIqIQrq9jlS47YeQits+Uvu1l5DWbtvW+HRve2N/3kd6OszQHH54+Pv5\nxaNtsxd1vnl0+ekjQkpFRCHV65HmzIZj2W5LPYR0vNxqTztYn4/a7a+HBDoh1e1cXkR9uDev\ny6yK82G/7HDJ5PqIp0Ph7ekM+7w57NferEr2kRBXSNc9m/zazenawW0ddHsd6XJaXrOG6m56\n5ecunu7d3o8fnIe4n4O06uw0tX9U5xPs2tKynxfZJqRURBXS6a9e2ou/86TvIdXrh2tITSJl\nW0FdQLHfNQt8N6S/y6bb072brF653M5scKv7TpF7Dul0rB9bXjYA11m9s/VthgkpFXGFFB1C\nSgUhSRFSKghJipBSQUhShJQKQpIipFQQkhQhpYKQpAgpFYQkRUipICQpQkoFIUkRUioISYqQ\nUhFTSJez3ornD6P7/Iag75dg7TXk+Me2jyekREQY0uVqDd2vf3jAj0uw9hpy/GPbxxNSIuIK\nqf3f+usbF948YEaElIoIQ+rfByHNOURaIg7pdk3Ul0ukXi63+nAJ1vru/PLdzXtsN52ptpds\nfZ7A+tTjsb9nmZAS4T2kf/28ndrDpt39mqjPl0i9vk325RKsl+/Obh8ncZ5qe8nWlwmsejy2\nzxNASEnwHVLPjt6WdDvY0FyWq3NN1KdLpHYvt3p+2OVaDIesudJCe4nW7f2KJZdLtnYmsLt8\ns/v52D5PACElIa410nnV0F7ernNN1KdLpJaPlxI63a8O1F4R6HzQr3P51fNBwIcJnL/Z/Xxs\nnyeAkJIQ3T5Snu0uf7ldTOHpyo7dRq5/vr14V/eb3kzgPumPj/09y4SUiOhC2p8vvkVI0xCS\nsehCqre0yvtfbl9/+mv3Vt+QnidASOgvvpAO54MNnWuiPl0itfiyj1R+DKkzgQ/7SG8e+3uW\nCSkR8YV0WSV1ron6dInU7uVWr5eSfDjy1pnW/VZnAt+O2j0+9vcsE1IiIgypOq+S7tdE7bzY\n014i9Xa51edLsN5fdXoN6WUC71+Denzs71kmpEREGNL1g8Rv10S9nX5wvUTq9XKr3UuwbrPb\n2QndaXVudSeQuWLfObPh02N/z/LiQrpcdRZPYgrp4yiaYXqfG/ttGgsLqXMBZzxYQEhH699r\n+wmzVfn902F7TmpZId1eI6CkZ/GHVG/HrWynePnM8y8fxNfbskK67ThS0ov4Qyrdm0/Um2Zb\nOJcbrI8WFtK9H0p6EX9IQVtaSJ2byf9uHxGS1JJCemiHkp4QktRiQ6KkJ4QktaCQnsshpAeE\nJLWskJ7/nvyvt4OQpJYT0ms2hNRFSFILDyn53+8dIUktOSRK6iIkqcWE9DYaQrojJKklhfT2\ni8n/hq8ISWrZIVHSHSFJLTwkNu5uCElqKSF9WvUQ0pU4pOQtJ6RPX6ekljYkLCOkz71Q0gUh\nqS0kpBF3JYWQ1AgpCYSktvyQKOlESHoLD4lV0hkhqS0gpK9rHUJqEZLaMkL6eiclEZLe0kNi\nldQiJLX4Q/qxzmGV1CAktUWENOn+JBCSGiElgZDUkgiJkghJbfkhUdKJkPQSCImNO0LSiz6k\nHusbQiIkuSWEZPAtS0dIaomElHpJhKSWQkiUREhysYfUrxFCIiSxBYRk9l0LRkhqhJQEQlJL\nJqS0SyIktTRCSn6VREhq6YSUdEmEpJZISKmvkghJjZCSQEhqkYfUe5Mt8W07QlKLPyTz71wi\nQlIjpCQQkhohJYGQ1FIKKeGSCEktmZDSXiURklpSIaVbEiGppRNS0qskQlIjpCQQklrcIQ3a\nXEt5246Q1KIPachUCQkqiYWUakmEpJZSSAmvkghJjZCSQEhqUYc0dFuNkKASe0jDJptsSYSk\nllRI6a6SCEktuZDSLImQ1NIKKdlVEiGpEVISCEmNkJJASGrphZRkSYSkllhIqa6SCEmNkJJA\nSGqElARCUvMZUrXO6j83uXPFn8UQY3Z4Et1JIiQ1jyEds3oZrjLXKgyGGNUEIUHBY0grV1b1\nH6tj3dTKracPMTakBEsiJDWPITlXXf6ot/JcNn2IcUkQEgS8hlT/kbnOXyYOQUi9EZKa1027\nw+m0af5o1khfd5IIyRghqXkM6eCy9eFUZnVJu9ztJg8xcm8nyZ0kQlLzefh7dzli19hMH2Js\nEIQEe35fkP1b5U1F5eZoMMSEkJIriZDUIj6zYXQPhARzhJQEQlIjpCQQktpcIRm8jjQlpNRK\nIiS1cEJyXf0mMf6lYUKCrRQ37RJcJRGSWpIhpbdKIiQ1QkoCIal5DWm/Kds9oHK9NxiCkPoj\nJDWPIVV552iCwRv7poWUVkmEpOYxpLXL/tpTv0/HXWbwxr4pMRASbHkMKTu/g6J1mP7Gvklr\nldRWSYSk5vuNfW//MmqIaSkQEkxFvUYaOSPTHx0dQlLzu4+0O799wmQfiZAGICQ1n4e/i85R\nu7yaOgQhDUBIan5fR1q3ryNl5cbgdaTJIaVUEiGpxXtmw8QSCAmWCCkJhKQWbUhTt80ICZZi\nDmnaEEmVREhqyYaU1iqJkNSSDimdkghJLd2QklolEZIaISWBkNQIKQmEpJZ2SMmUREhqCYeU\n0iqJkNQSDymVkghJLeWQElolEZIaISWBkNQIKQmEpJZ6SImUREhqSYeUziqJkNSSDymNkghJ\nLdaQjBIgJNiIOCSTgQgJJgjJYjLBIyQ1QkqiJEJSSzykVFZJhKRGSEmUREhqqYeUyCqJkNQI\niZBggJAICQYIiZBgINKQ7I4RpHG0gZDU4g3JbCxCwnSEREgwQEiEBAOElMROEiGpEVISqyRC\nUiOkJFZJhKRGSEmskghJjZAICQYIiZBggJCS2EkiJDVCMp9aiAhJjZBOKaySCEmNkASTCw8h\nqRGSYHLhISQ1QhJMLjyEpEZIl8ktuyRCUiMkyfRCQ0hqhHSd3qJLIiQ1QhJNMCyEpBZnSPYr\nEELCJNGGZD0iIWEKQrpNccklEZIaIekmGRBCUiMk3SQDQkhqhKSbZEAISY2Q7pNccEmEpEZI\nymkGg5DUCKkzzeWWREhqhKSdaCAISY2QtBMNBCGpEZJ2ooEgJDVC0k40EISkRkjdiS62JEJS\nIyT1VINASGqE9DDVpZZESGqEpJ9sAAhJjZD0kw0AIakRkn6yASAkNUJ6nOxCSyIkNULyMd3Z\nEZIaIT1Nd5klEZIaIfmZ8MwISY2Q/Ex4ZoSkFmVIui0wQsI4sYakGnihJRGSGiF5m/ScCEmN\nkF4mvcSSCEmNkDxOez6EpEZIHqc9H0JSIySP054PIakR0uu0F1gSIakRkteJz4WQ1AjJ68Tn\nQkhqhPRm4ssriZDUCMnz1OdBSGqE9G7qiyuJkNQIyfvk50BIajGGpF5lEBIGizQk7eiEhKEI\n6e30l1YSIakR0hwDeEdIaoT0foCFlURIaoQ0zwieEZIaIc0zgmeEpEZI84zgGSGpEdKHEZZV\nEiGpEdJcQ3hFSGqE9GmIRZVESGoTQ8o3R7NZ+TDEm3s9LOSEhEEmhlQvb4qWCMkYIalNDKn6\nWylaIiRjhKRmsI+03+TWLQUR0pJKIiQ1m4MNh6xe7rbT5+bLEN17fSzjhIQhTELaFc0/4K4w\nmJ9PQzzcS0hDEZLa9JCqTb06yndVXVNpM0+EZI6Q1KaGtG8ONqwP5zvMFr0wQlpQSYSkNvV1\npHpltK2ud2QWc/Q8xJt7vSzihIQBpr6OVO7MZuXDEG/u9RXSckoiJLWpryOZzcjHId7c62cJ\nJyT0N/nMhsuNzGyz7nmIN/cS0lCEpGYU0tF2oQslpMWUREhqE0Laua7c31z5WsAJCb1NWSPl\n3Y72/ubKY0hLKYmQ1Kz2kWyFEdKCVkmEpMYb+0IYSI6Q1CaE1CxlnY07f3NFSIMRkhohfR9o\nISURkhqbdmGMJEZIarOE9HPx/Hq/x/XEYlZJhKQ2NaRtfjod84FHvyeHNGSwSQgJ/UwMadcs\naM3bY93vktyj0XNFSMMRktrEkAr3dzq4/PTX4+2x+yzKkBZREiGpGbwge3Drfq/MVqUrjtdH\njZ8rr8s2IaEXg5BKt+t7isOfq9dgsYW0hJIISW3ypt1h17wxts+mXeNYuLKKKaSFrJIISW36\nwQbnNs3i1vedshuX7QjJN0JSm3z4O2v2kE75X++HH/LfyyYhGSMktTlekF1FFtICSiIktQhP\nEfK8ZBMSeiCkXwgJPUwNaZP3eYH1zUTieEH2tJBtO0JSmxjSZuzbKF6/v/9pD54XbELCbxND\nyiw/g+L9EK93+g8p+pIISS3CazZ4X64JCT9NDKl0kmutEpIxQlKbGNIxKwa8E2m/Kds9oHL9\n40GEZIyQ1CZv2vU/2FB1r4P3/dS8bxPzv8uygJ0kQlLzGNLaZX/nD1I67s4nFo2aqxmWakLC\nLx5fkM3c4Xb78P2zlMILKfKSCEnNY0jOffrLoCHmWKgJCT9MDmlXtm/uO/5+XLxrJELCL1ND\nKs6bPS77XVK9j7Q7f1ds+0jxb9sRktrEkLauqJplbOtWvx9YdA5N5F9ffwotpOhXSYSkNvkU\nocsbx3staPt1+zpSVm4mvI40V0hRl0RIaganCPUPadQQL/fNskgTEr6aGFJ+WSMd/H1iHyGN\nQEhqNvtIO+OzwIMMKeaSCElt6lG7stcpP5OGeL5vniWakPCNyetIrux/EaERQzzdN1tIEZdE\nSGrxXbNhrgWakPAFIfVFSPhiWki7VfPWiOLX24umDPF633whxVsSIalNCel4P1Wh6HGu3Zgh\n3t031/JMSPhsQkhV5vJdc6bP8S//fg7q6CHe3kdIwxGS2oSQ1p1j3kVzJX07hGSMkNQmhJS7\n+/bc0faFpFBDirYkQlKbENKAN+qNHeLtfbMtzhGXREhqhDQAIeETQhqAkPAJIQ1ASPhkUkgP\nfM3VvCFFWhIhqRHSENGWREhqnGs3CCHhPUIahJDwHiENEuu2HSGpEdIwhIS3CGmYSFdJhKRG\nSAMREt4hpIEICe8Q0kBxbtsRkhohDTX3+KMQkhohDTX3+KMQkhohDRXlth0hqRHSYLPPwAiE\npEZIg8W4SiIkNUIabv45GIyQ1AhpuPnnYDBCUiOk4SLctiMkNUIaIYBZGIiQ1KILKYTVQQjz\nMAwhqcUYkn78X4KYiSEISY2QxghiJoYgJDVCGiO6bTtCUiOkUWIriZDUCGmcQGajL0JSI6Rx\nApmNvghJjZDGiWzbjpDUCGmkUOajH0JSI6SR4lolEZIaIY0VzIz0QUhqhDRWMDPSByGpEdJY\nUW3bEZIaIY0Wzpz8RkhqhDRaOHPyGyGpEdJoMW3bEZIaIY0XUUmEpEZIE4Q0L98RkhohTRDS\nvHxHSGqENEE823aEpEZIU0RTEiGpEdIkYc3NZ4SkRkiThDU3nxGSGiFNEsu2HSGpEdI0gc3O\nJ4SkRkjTRLJKIiS12EIKbsENbX7eIyS1CEPSDz9EcGW/RUhqhDRVeHP0BiGpEdJU4c3RG4Sk\nRkhTRbFtR0hqhDRZgLP0gpDUCGmyGFZJhKRGSNOFOE9PCEmNkKaLYJVESGqEZCDImXpASGqE\nZCD8VRIhqRGShTDnqoOQ1AjJQvCrJEJSIyQTgc7WDSGpEZKJ0FdJhKRGSDZCna8LQlIjJBuB\nr5IISY2QjAQ7Yy1CUiMkI2GvkghJjZCshDtnJ0LSIyQrQa+SCEmNkMwEPGuEJEdIZgKeNUKS\nIyQzIW/bEZIaIdkJuCRCUiMkQ+HOHCGpEZKhcFdJhKRGSJaCnTtCUiMkS8GukghJjZBMhTp7\nhKRGSKZCXSURkhoh2Qq0JEJSIyRjYc4gIakRkrEwZ5CQ1AjJWJjbdoSkRkjWgiyJkNQIyVyI\ns0hIaoRkLsRVEiGpEZK9AOeRkNQIyV6A80hIaoRkL8BtO0JSIySB8EoiJDVCUghuLglJjZAU\nglslEZIaIUmEVhIhqRGSRmDzSUhqhKQR2CqJkNQISSSsGSUkNUISCWuVREhqhKQS1JwSkhoh\nqQS1SiIkNUKSCWlWCUmNkGRCWiURkhoh6QRUEiGpEZJQODNLSGqEJBTOKomQ1AhJKZiSCEnN\nZ0jVyrlid5nI16ksJaRgZpeQ1DyGVGXNP9CuPE8kmZCCmF9CUvMY0tpt65q2WdFOJI2QQimJ\nkNQ8hpSdH3jM8mM6IQUyw4Sk5jGk6xJVFUVaIQUwx4Sk5jGk3FXXW0U6IYVREiGpeQxp61aX\nW0dXpBNSELNMSGo+D3+vb0vU7sfCtbSQZp9nQlLz+oLsobzeOq7SCSmEeSYkNc5s0AtglURI\naoTkwfwlEZIaIfkw+1wTktpcISV0sOEUwCqJkNTCCcl1fXlcjCHNXhIhqbFp5wchLVxkIc39\nL/toM884IanFF5J+dIl5SyIkNa8h7Tfl+S1J6/3IIeINad5ZJyQ1n2/syztHE4pxQ8Qd0nzz\nTkhqXt/Yl/0d2lvHXebWo4aIOKRZ552Q1Ly+se9wu31w2aghIg9ptpknJLUZ3tj3+pf+Q8Qc\n0pwlEZIaaySP5iuJkNT87iPtju2tNPeRLidvzDIyIan5PPxddI7a5dW371xoSPOtkghJze/r\nSOv2daSs3CT4OlJrrpIISY0zG/yaqSRCUiMkzwhpmQjJs3lWSYSkRki+zVISIakRkndzlERI\naoTkHyEtECH5N8MqiZDUCGkG/ksiJDVCmoP3kghJjZDm4P2kO0JSI6RZENLSENI8PJdESGqE\nNBO/JRGSGiHNxWtJhKRGSLPxWRIhqRHSfAhpQQhpPh5XSYSkRkgz8lcSIakR0py8lURIaoQ0\nK18lEZIaIc3LU0mEpEZIM/NTEiGpEdLM/Jy/SkhqhDQ3QloEQpqdj5IISY2Q5uehJEJSI6QA\n6EsiJDVCCoG8JEJSI6QgqEsiJDVCCoO4JEJSI6RAaEsiJDVCCoW0JEJSI6RgKEsiJDVCCoew\nJEJSI6SA6EoiJDVCCojuBFZCUiOkkMhKIiQ1QgqKqiRCUiOksIhKIiQ1QgqNpCRCUiOk4ChK\nIiQ1QgqPoCRCUiOk8Aj2kwhJjZACZF8SIakRUpCsSyIkNUIKk3FJhKRGSIGy3bwjJDVCCpVp\nSYSkRkjBsiyJkNQIKWB2JRGSGiGFzKwkQlIjpKBZlURIaoQUNqMdJUJSI6TA2ZRESGqEFDqT\nkghJjZDCZ1ASIakRUgSmr5QISY2QYjC5JEJSI6QoTC2JkNQIKRLTUiIkNUKKxaSSCEmNkKLh\nJqRESGqEFJHxJRGSGiHFZHRJhKRGSHEZmRIhqRFSZMbtKRGSGiHFZlRJhKRGSPEZkRIhqRFS\njAaXREhqhBSloSURkhohRWpYSoSkRkixGrSnREhqhBStISURkhohRax/SoSkRkhR65sSIakR\nUtxcv5QISY2QYterJEJSI6T49UiJkNQIaQl+pkRIaoS0CL92lQhJjZAW4ntKhKRGSIvxLSVC\nUiOkBfmcEiGpEdKiuA8tEZIaIS3M+5IISY2QFuddSoSkRkhL9NISIakR0iI97ysRkhohLdRj\nSoSkRkjL1WmJkNQIacmuKTlCUiOkZSMkTwhp8S4h8cxJEVIKziHx3AkRUhIuIfH0yRBSEpp9\npNuRh7lnZpEIKQnngw2OlmQIKQm3o3a0JEJISege/qYlBUJKwtPrSLRkjpCS8PqCLC3ZIqQk\nvD2zwRGTHUJKwqdThGjJCiEl4cu5do6YLBBSEn6ctEpMkxFSEn6f/e0cNU1BSEno9zYKWhqP\nkJLQ//1IrJnGIaQkDHxjHzUNRkhJGP4OWeeoaQhCSsLIt5pTU2+ElIQp12ygpj4IKQmTL37y\nXBO/hieElASbqwi91MRv44aQkmB4Oa43NfE7IaREmF/X7l1OSf9qCCkJqgtEvu0pyV+R15D2\nm7J9osv1fuQQif6WJtNeafV9Tmn9qjyGVOWd57gYN0Rivx0zXi5Z/KmnJH5pHkNau+zv0N46\n7jK3HjVEEr8TAa/X/v4S1HJ/fR5DytzhdvvgslFDLPg3ITXLRfS/BrW036THkB6eue9PIyHZ\n+lebc/xfSS3g18oaafn+Xcw9H40eScWZld99pN2xvcU+kkf/Ouaelwc9m4okLZ+Hv4vO85JX\no4aI4BkNzL9/wZbUNayqAOPy+zrSun0dKSs3vI7kSyQhPRqV1byJRXZmg4/BF+XfvyhLejE5\nLXVrsYWEYZYS0jv2cU2oj5AW7bmjZZX00QwlzRXSyNeRMMxrSBAZ/ssRhdRzFYoh5l66EjL8\nl8OmXTzmXroSMvyXQ0gRmf7rRi8cbFg2QvIksjf2YSg68iOyN/ZhKELyI7I39mEwOvIisrdR\nYAQy8iCyN/ZhFDKSY42UhFneap6UyN7Yh3EISS2yN/ZhHEJSi+yNfRiHkNQ4syEJhKRGSEkg\nJDVCSgIhqRFSEghJjZCSQEhqhJQEQlIjpCQQkhohJYGQ1AgpCYSkRkhJICS1QEOCrTqkuWdh\n6UYs5fbhhDj2Ioda5A8V6/NHSPEOtcgfKtbnj5DiHWqRP1Sszx8hxTvUIn+oWJ8/Qop3qEX+\nULE+f4QU71CL/KFiff4IKd6hFvlDxfr8EVK8Qy3yh4r1+SOkeIda5A8V6/NHSPEOtcgfKtbn\nj5DiHWqRP1Sszx8hxTvUIn+oWJ8/Ts0GDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBAS\nYICQAAOEBBggJMAAIQEGCAkwQEiAgdlCWmcuW1eeBtv6+jG3uaefqlo5tzp4GOhs7+cJHH0J\n++EOzRN4tJveXCEV7TOW+xns4Od3U//r0P5UmYeSsnYkXyVVmZ+F219IO+vf1Ewh7V12OB0y\nt/cxWD2Onx/z4FZVs/pbyUdaN2OsXSkf6Kz0tJbw9gOdsnr5q0q3NpvgTCGt3a7+889tPIy1\ndYWnkMrzMB5Gy1zlZ6DWn6fNra2X5aHx1yZUucxsijOFVLpm89TPv0D1c+ZribsM6Gs0w+Xg\nm6Ovf4m2butjmNrKfKt4ppCct3+761w9LtqNyhV+Blp7Wu4Kd/TzBJZut3KZ3ebWZ7k7bbJ2\nQ9xKAiF5HKe1bTdb5ertLR+L3Om0cX+ensDyfKzBw79DzrVjGa7RCcnaMfOzx7wtMy+7FO32\nt58n0NXJniofK1rXHOyqVoZPICEZqzJPG3anZkvfw7Zd3hwj9vkvUeXhVZHzKwdHw5FmCilb\nbEiFp/vpHMwAAAPcSURBVNfGGpZHnT5ZtRuqXreNPQxm/w/5rEftjr5eN/C2HBzzwvDV8p+8\nLHLO4/kGlyHlQ9i/UDFTSJv237mdp91lbyHtfB2wO7+OZLlp8onPkK4/lf6f1/PydzT8daVw\nZoO3kCx/Md+1ZzZUpbfXXfw8gevmH9Zq7eGoZ/1PUNUcbPgzm+Jc59rlvg50tjyFtPL5j7fP\n58/TE1idfyofmykb6+dvrpCq9uxvX6N5Csnn7kT9/OW+1kfensDK30+1K2yXP96PBBggJMAA\nIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAA\nIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggpLns\nO5+knD98bOqAj8erVvdPivz6eYEvH8uq/5zWtBDSbO4fRX10WfeOASGVdTeb28M+h5Q/f/Hl\nC5iG53M2m9u6ZP34+cMDQnLu2OthL/d5+lDYdPB8zqa6rYeyTg6ngSH1exghqfF8zqe47Bnt\n24+p39WbaefP2W4W8vOCfv5zm7vs4aO+6y+0n/39uCX3po1d4Vyxu3xjdwxfn72eDp7N+ezc\nqv3/qglqc67ivJQ/hFS2dxT3xxXXL/wKaXv+hu21m/sYhGSNZ3NG2fnZbxdp5/5Op7/LzW5I\nO1dUp6q4H9f7c9nhdMia73/atHs+1pC5Q/Pt+X1yD2PAEM/njNZtDH+dQw1vQipddWp2qMrr\n95RtU7t2HfU9JHer72VfipCM8XzO6NDGUDTrjdpxtynehPRyVNs9rMfuU3ttY+1ceTh073sY\nA4Z4PueU12ub6vLCbHHLZXJIt+/fZPX/s+PtvscxYIjnc05btzltXHtEbuXy7e74NqSnBw0I\nqd4CXOf3faSnMWCI53NOzdoob/eBzkv2U0j78z7S09k8132k8vRr0+7+9c5BQEKS4Pmc1cpd\nj4E7tz8d7vsveb2eqtq/tgfp6nXX7WDD56N2z1PPz0fp8usZEA9jHJ+/G1MQ0qx27npkbX3Z\nE9qfi2hfAirbZf68Y5PdF/zb60i/Qvq7TTJ3zWkUnTHOX4AZQppXdlueV3Ub+3aD7XqcYHU/\ns8GtuiuQbXY+s+Hnpl17ZkNzauw+b8e5j3H5AqwQEmCAkAADhAQYICTAACEBBggJMEBIgAFC\nAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFC\nAgwQEmCAkAADhAQYICTAACEBBggJMEBIgIH/AeRKEpqeDME6AAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualization of critical value, rejection region for F-test\n",
    "r=4;n=33;alpha=0.05;\n",
    "x.grid=seq(from=1e-5,to=6,length.out=1000);\n",
    "density.grid=df(x=x.grid, df1=r-1, df2=n-r)\n",
    "critical.value=qf(1-alpha,df1=r-1,df2=n-r);\n",
    "plot(density.grid~x.grid,type='l',xlab=\"Value of F-stat\",ylab=\"Density\",lwd=3,xlim=c(0,6),ylim=c(0,0.8))\n",
    "abline(v=critical.value,lwd=3,col='red')\n",
    "segments(x0=critical.value,x1=10,y0=0,y1=0,lwd=3,col=\"orange\")\n",
    "points(x=critical.value,y=0,pch=16,col=\"orange\",cex=2)\n",
    "legend(x=3.8,y=0.8,legend=c(paste0('Critical value ', round(critical.value,digits=2)), 'Rejection region'),lty=1,lwd=3,col=c('Red','Orange'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Example.** In the Spock trial, we want to test the null hypothesis $H_0: \\mu_1=\\mu_2=\\mu_3=\\mu_4$ against the alternative $H_1:$ not all $\\mu_i$'s are equal. We can calculate the F-statistics $F^*={\\rm MSTR}/{\\rm MSE}=17.6$, when $F(0.95;3,29)=2.93$. We can thus reject the null hypothesis at the nominal significance level $0.05$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Rationale.** In real life, we reject a hypothesis if the null is *unlikely to be true* given our observations. We say the null hypothesis is *unlikely to be true* if the observed data is very **extreme** under the null hypothesis. We often calculate the p-value, which is *the probability of observing a more extreme test statistic* $F$ *than the test statistic* $F^*$ *calculated on the current sample.* In other words,  \n",
    "$$ {\\rm p value} \\ =\\mathbb{P}(F \\geq F^* ).$$\n",
    "We should see the p-value as a summary measure of the plausibility of the null hypothesis. Noting the subjectivity in \"plausibility\", we thus should specify the significance level (i.e., the acceptable type I error rate) before data analysis. \n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Distribution of $F^*$ under the alternative.** The numerator ${\\rm SSTR}/\\sigma^2$ has a non-central Chisquare distribution with d.f. $r-1$ and non-centrality parameter $\\sum_{i=1}^r n_i (\\mu_i-\\mu)^2/\\sigma^2$. We know that $F^*$ follows a non-central $F$-distribution with $df=(r-1,n_T-r)$ and the noncentrality parameter $\\sum_{i=1}^r n_i (\\mu_i-\\mu)^2/\\sigma^2$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "**Expectation of $F^*$.** To derive the expectation of $F^*$, we have \n",
    "$$\n",
    "\\mathbb{E}[F^*]=\\mathbb{E}\\left[\\frac{\\rm MSTR}{\\rm MSE} \\right]=\\mathbb{E}\\big[{\\rm MSTR}\\big]\\mathbb{E}\\left[\\frac{1}{\\rm MSE} \\right], \n",
    "$$\n",
    "where the second equality holds due to independence between MSTR and MSE. We know that the numerator $\\mathbb{E}[{\\rm MSTR}]=(n_T-r)^{-1}\\sum_{i=1}^r n_i (\\mu_i-\\mu)^2/\\sigma^2$. It remains to derive the expectation of the denominator. \n",
    "\n",
    "We have \n",
    "$$\n",
    "\\left[ \\frac{1}{\\rm MSE} \\right]=\\frac{1}{\\sigma^2} (n_T-r) \\mathbb{E}\\left[\\frac{1}{\\chi^2_{n_T-r}} \\right]=\\frac{1}{\\sigma^2} \\frac{n_T-r}{n_T-r-2},\n",
    "$$\n",
    "because $(n_T-r) {\\rm MSE}/\\sigma^2 \\sim \\chi^2_{n_T-r}$ and $\\mathbb{E}[1/\\chi^2_{v}]=1/(v-2)$ for $v>2$. \n",
    "\n",
    "Therefore,\n",
    "\\begin{align}\n",
    "\\mathbb{E}[F]= &  \\frac{n_T-r}{n_T-r-2}\\left[ 1+\\sum_{i=1}^r \\frac{n_i}{r-1} \\left(\\frac{\\mu_i-\\mu}{\\sigma}\\right)^2 \\right]\\\\\n",
    "\\approx &  1+\\sum_{i=1}^r \\frac{n_i}{r-1} \\left(\\frac{\\mu_i-\\mu}{\\sigma}\\right)^2.\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "**Interpretation of $F^*$.**\n",
    "1. ${\\rm MSE}$ is an estimator of $\\sigma^2$, but the maximum likelihood estimator for $\\sigma^2$ is ${\\rm SSE}/n_T$.\n",
    "3. The ratio $F^*$ fluctuates about \n",
    "$$\n",
    "\\frac{\\mathbb{E}[{\\rm MSTR}]}{\\mathbb{E}[{\\rm MSE}]}=1+\\sum_{i=1}^r \\frac{n_i}{r-1} \\left(\\frac{\\mu_i-\\mu}{\\sigma} \\right)^2,\n",
    "$$\n",
    "where the second term is a unit-free measure of the variablity among $\\mu_1,\\mu_2,\\ldots, \\mu_r$. Under the null, the second term equals zero and thus $F$ fluctuates around one. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 4.1.7 Testing linear combination \n",
    "\n",
    "For any $c_i, i=1,\\ldots, r$, the least squares and maximum likelihood estimator of $L=\\sum_{i=1}^r c_i \\mu_i$ is $\\hat{L}=\\sum_{i=1}^r c_i \\bar{Y}_{i\\cdot}$, and further \n",
    "$$\n",
    "\\mathbb{E}\\big[\\hat{L} \\big]= L, \\ {\\rm and}\\  {\\rm var}\\big(\\hat{L}\\big)=\\sum_{i=1}^r \\frac{c_i^2}{n_i}\\sigma^2.\n",
    "$$\n",
    "Therefore, an unbiased estimator of ${\\rm var}\\big(\\hat{L}\\big)$ is \n",
    "$$\n",
    "s^2\\big(\\hat{L}\\big)=\\sum_{i=1}^r \\frac{c_i^2}{n_i}{\\rm MSE}. \n",
    "$$\n",
    "Under the null hypothesis $H_0: L=0$ and the normality assumption on $\\epsilon$, we have  \n",
    "$$\n",
    "\\frac{\\hat{L}-L}{ s\\big(\\hat{L}\\big)} \\sim t(n_T-r).\n",
    "$$\n",
    "\n",
    "To test the hypothesis $H_0: L=0$ against $H_1: L \\neq 0$. We can calculate the t-statistic $t^*=\\hat{L}/s(\\hat{L})$. We can calculate the quantile of the t-distribution as before to finish the test.\n",
    "\n",
    "**Example.** Consider the quantity $L=\\mu_1-(\\mu_2+\\mu_3+\\mu_4)/3$. We can calculate that $\\hat{L}=-15.217$, and  $s^2(\\hat{L})=4.648$. We can then test the null hypothesis or construct a $99\\%$ confidence interval with $t(1-0.01/2;n_T-r)=2.756$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "15.217"
      ],
      "text/latex": [
       "15.217"
      ],
      "text/markdown": [
       "15.217"
      ],
      "text/plain": [
       "[1] 15.217"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "4.648"
      ],
      "text/latex": [
       "4.648"
      ],
      "text/markdown": [
       "4.648"
      ],
      "text/plain": [
       "[1] 4.648"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "2.756"
      ],
      "text/latex": [
       "2.756"
      ],
      "text/markdown": [
       "2.756"
      ],
      "text/plain": [
       "[1] 2.756"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "(Lhat=unname(round(sum(anova.fit[[1]][2:4])/3,digits=3)))\n",
    "(s2=round(sum(anova.fit[[2]]^2)/anova.fit[[8]]*sum(c(1,-1/3,-1/3,-1/3)^2/table(Spock[,1])),digits=3))\n",
    "(tquantile=round(qt(1-0.01/2,anova.fit[[8]]),digits=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.8 Simultaneous inference[就是牺牲conf.interval最紧性质，移动C.I.两端，使得对不同L，这个新的C.I.都是“95%”的]\n",
    "\n",
    "Assume that we have a few linear combinations, $L_1, L_2, \\ldots, L_m$. Let $L\\equiv (L_1,\\ldots, L_m)^T$. We can construct a simultaneous interval for $L$ as $\\hat{L}\\mp K s\\big(\\hat{L}\\big)$ for some appropriate multiplier $K$\n",
    "\n",
    "[eg:$K:=$(1-$\\alpha/2$) quantile of the distribution, remember it is a two-sided one], \n",
    "\n",
    "where $\\hat{L}$ is an estimator of $L$ and $s(\\hat{L})$ is the standard deviation of $\\hat{L}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mse=sum(anova.fit$residuals^2)/anova.fit$df.residual;\n",
    "\n",
    "# Create vectors for the first two linear combinations \n",
    "comb.mat<-matrix(0,nrow=2,ncol=4)\n",
    "comb.mat[1,]=c(1,-1,0,0);comb.mat[2,]=c(1,0,-1,0);\n",
    "\n",
    "# Obtain the estimates \n",
    "diff = numeric(dim(comb.mat)[1]);\n",
    "diff.sd=diff;\n",
    "mean.tmp=anova.fit$coefficients;mean.tmp[1]=0;\n",
    "ns=as.numeric(table(Spock$Judge));\n",
    "for(i in 1:length(diff)){\n",
    "  diff[i]=sum(comb.mat[i,]*mean.tmp);\n",
    "  diff.sd[i]=sqrt(sum(comb.mat[i,]^2/ns)*mse);\n",
    "}\n",
    "alpha=0.05;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Bonferroni method【缺点是显著度会不断下降】. The multiplier is $K=t\\big(1-\\alpha/(2m); n_T-r\\big)$. Then, the simultaneous intervals cover the true parameters with probability at least $(1-\\alpha)$.\n",
    "\n",
    "$P(A_1^C)=0.05,P(A_2^C)=0.05$.\n",
    "\n",
    "$P(A_1^C\\cup A_2^C)\\le P(A_1^C)+P(A_2^C)=0.1$.\n",
    "\n",
    "$P((A_1^C\\cup A_2^C)^C)\\ge 0.9$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.83155268751868"
      ],
      "text/latex": [
       "2.83155268751868"
      ],
      "text/markdown": [
       "2.83155268751868"
      ],
      "text/plain": [
       "[1] 2.831553"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Bonferroni correction:\n",
    "m=6; # for all pairwise differences, although we only show two here\n",
    "B.stat=qt(1-alpha/(2*m),anova.fit$df.residual);\n",
    "B.stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Tukey-Kramer method【仅仅在做pairwise $\\mu_i -\\mu_{i'}$ 比较时有效，在其他情况下如$\\mu_i -0.5\\mu_{i'}-0.5\\mu_{i''}$都无效】. This approach only works for pairwise comparisons, e.g., $\\mu_i -\\mu_{i'}$. The $100(1-\\alpha)\\%$ confidence interval for $\\{\\mu_i-\\mu_{i'}: i, i' \\in \\{1,\\ldots, r\\}, i\\neq i'\\}$ is \n",
    "$$\n",
    "[\\bar{Y}_{i\\cdot} - \\bar{Y}_{i' \\cdot}] \\mp T s\\big([ \\bar{Y}_{i\\cdot} -\\bar{Y}_{i'\\cdot}] \\big), \\ i\\neq i', \\ T=\\frac{1}{\\sqrt{2}} q(1-\\alpha; r, n_T-r), \n",
    "$$\n",
    "where $q$ is the studentized range distribution. The coverage is exactly $1-\\alpha$ for a balanced ANOVA model, and at least $1-\\alpha$ for unbalanced cases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".dl-inline {width: auto; margin:0; padding: 0}\n",
       ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
       ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
       ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
       "</style><dl class=dl-inline><dt>(Intercept)</dt><dd>14.6222222222222</dd><dt>Judge2</dt><dd>18.9944444444444</dd><dt>Judge3</dt><dd>14.4777777777778</dd><dt>Judge4</dt><dd>12.1777777777778</dd></dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[(Intercept)] 14.6222222222222\n",
       "\\item[Judge2] 18.9944444444444\n",
       "\\item[Judge3] 14.4777777777778\n",
       "\\item[Judge4] 12.1777777777778\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "(Intercept)\n",
       ":   14.6222222222222Judge2\n",
       ":   18.9944444444444Judge3\n",
       ":   14.4777777777778Judge4\n",
       ":   12.1777777777778\n",
       "\n"
      ],
      "text/plain": [
       "(Intercept)      Judge2      Judge3      Judge4 \n",
       "   14.62222    18.99444    14.47778    12.17778 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "29"
      ],
      "text/latex": [
       "29"
      ],
      "text/markdown": [
       "29"
      ],
      "text/plain": [
       "[1] 29"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "4"
      ],
      "text/latex": [
       "4"
      ],
      "text/markdown": [
       "4"
      ],
      "text/plain": [
       "[1] 4"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "2.72450905738992"
      ],
      "text/latex": [
       "2.72450905738992"
      ],
      "text/markdown": [
       "2.72450905738992"
      ],
      "text/plain": [
       "[1] 2.724509"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Tukey-Kramer \n",
    "T.stat=qtukey(1-alpha, nmeans=length(anova.fit$coefficients), df=anova.fit$df.residual)/sqrt(2);\n",
    "anova.fit$coefficients\n",
    "anova.fit$df.residual\n",
    "length(anova.fit$coefficients)\n",
    "T.stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Scheffe method【对所有有限的c都有效,$L_j=\\sum_{i=1}^r c_{ji} \\mu_i$ and $\\sum_{i=1}^r c_{ji} =0$】. This approach applies to all possible constrasts, where the coefficients sum up to zero for each $L_j$, i.e.,$L_j=\\sum_{i=1}^r c_{ji} \\mu_i$ and $\\sum_{i=1}^r c_{ji} =0$. The multiplier is $S= \\big[(r-1) F(1-\\alpha; r-1,n_T-r)\\big]^{1/2}$. The coverage is exactly $1-\\alpha$ for _all_ possible contrasts. For finitely many contrasts, the coverage is at least $1-\\alpha$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.96683158756821"
      ],
      "text/latex": [
       "2.96683158756821"
      ],
      "text/markdown": [
       "2.96683158756821"
      ],
      "text/plain": [
       "[1] 2.966832"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Scheffe \n",
    "S.stat=sqrt( (length(anova.fit$coefficients)-1)*qf(1-alpha,length(anova.fit$coefficients)-1,anova.fit$df.residual))\n",
    "S.stat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Choosing which method**\n",
    "- All pair comparison, (near) balanced design, then Tukey method is the best.\n",
    "- If testing $m$ hypothesis, $m$ is small use Bonferroni; If $m$ large, use Scheffe method. \n",
    "- 在所有方法都合法时，选一个CI最紧的！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A matrix: 1 × 3 of type dbl</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>Bonferroni</th><th scope=col>Tukey</th><th scope=col>Scheffe</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>2.831553</td><td>2.724509</td><td>2.966832</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A matrix: 1 × 3 of type dbl\n",
       "\\begin{tabular}{lll}\n",
       " Bonferroni & Tukey & Scheffe\\\\\n",
       "\\hline\n",
       "\t 2.831553 & 2.724509 & 2.966832\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A matrix: 1 × 3 of type dbl\n",
       "\n",
       "| Bonferroni | Tukey | Scheffe |\n",
       "|---|---|---|\n",
       "| 2.831553 | 2.724509 | 2.966832 |\n",
       "\n"
      ],
      "text/plain": [
       "     Bonferroni Tukey    Scheffe \n",
       "[1,] 2.831553   2.724509 2.966832"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table.stats=matrix(0,1,3);\n",
    "table.stats[1,]=c(B.stat,T.stat,S.stat);\n",
    "colnames(table.stats)=c('Bonferroni', 'Tukey', 'Scheffe')\n",
    "table.stats\n",
    "# Then, we can construct the confidence intervals as, e.g.,\n",
    "CI.bonferroni =matrix(0,nrow=2,ncol=2);\n",
    "for(i in 1:length(diff)){\n",
    "  CI.bonferroni[i,]=diff[i]+c(1,-1)*B.stat*diff.sd[i];\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.9 Diagnostics of one-way ANOVA\n",
    "\n",
    "Recall that the (cell means) model for a one-way ANOVA is \n",
    "$$\n",
    "Y_{ij}=\\mu_i +\\epsilon_{ij}, \\  j=1,\\ldots, n_i, \\ i=1,\\ldots,r,\n",
    "$$\n",
    "where $\\epsilon_{ij}$ are i.i.d. $N(0,\\sigma^2)$. \n",
    "\n",
    "**Possible departures** from the model assumptions\n",
    "1. Variances of error terms are unequal.\n",
    "2. Error terms are not independent.\n",
    "3. Error terms are not normally distributed.\n",
    "4. Outliers (samples that do not follow this model).\n",
    "5. Missing variables. \n",
    "\n",
    "All diagnostics rely on the _residuals_.\n",
    "- Residuals $e_{ij}\\equiv Y_{ij}-\\bar{Y}_{i\\cdot}$.\n",
    "- Semistudentized residuals $e^*_{ij}=e_{ij}/\\sqrt{\\rm MSE}$\n",
    "- Studentized residuals $r_{ij}=e_{ij}/s(e_{ij})$, where $s^2(e_{ij})=(1-1/n_i) {\\rm MSE}$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDT09PZ2dnh4eHp6enw8PD///8uNL8wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAfe0lEQVR4nO3d20LiWBBG4R1OIgOE93/aIeEYQKtN/iR7V63vYtpuW4ok\ntUZFWtMJwGBp7jsAeEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBI\ngAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAA\nIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASP2llLovPf7g2XqSO/NVpfSn\nSS939vN9/+0t8IxT098/hbSvJjnFX+fZhDQjTk1//xTSRNu3SOnwt7cgJClOTX9vIf3+lya6\nL+PdACH9glPT30/vkeqv5fml1fflT9Ltr+3WzUdfu+ubHM+/W26f3vK4SJvzS9+r88uLzfF2\ne9tFWuxPp22Vlvvu+M7tvb9/ud5evalSdbm55zv2eIvjukqL7es71duvr3fn9VZwRUj9/RDS\nsbrms+yEtLy+vGrfYn/9K4+3XLRvcPtbaX/508vvj5v7n909395zr6fO7d3uzP7ljt3v8/sd\neT6i17vzdiu4IqT+fgjp/I7i/H/r+ryD26cNX91W8lJSdf/t7S1T82bb837Wp9Om00dK1XOD\nF53b+xRSe3v3OdXLHbvf5/c78nREb3fn7VZwRUj9pWfXP7j8t/lIqD6/U7j/2Wl3/nVbnz8s\nOv96/mjs+7zazS/V4y2bjW0eMzh2bun8p9vm3cuh/eUx++X2Pnxo197epYT6svyf7tj1juyq\njyG93Z23w8MVIfX3Q0jNTt4/Fbpt5vr2P/BN+yj16rL+zR7f3nL3ctOX/+47vzz+wsvtfQip\nvb3znPry+9XnO3a7I7uPIb3dnbfDwxUh9fdDSF+XP7gu2+NV7Uafju0fVLc9fX31+S98b5bp\neXM/Prj+cnsfQqq7d7H6fMfS+x15/u2Hu9O9FVwRUn+P5e1u4ea2u8e3V91e+ml/T9+LpzJ/\nD6nz0ntIt1+fbu/DHTNC+nR3OreCK0Lq76eQTvX35dGuZedV9/cg1cf3SO1vmw/1Fuvt4U/v\nkarXVz5+W3X/+P2O/R7Sx7vTuRVcEVJ/P4bUaL/K8/izlfk5UvvaxfXPzZBW5udIt7/2+iFY\n9459uCNNoPt0e7Dh/e50bgVXnIz+fghpcf/85Pauov7xUbv0+SMq+z2S/ahd+2szZ9/+svx0\nx9qH9ToPH55/3bRPEHy+A527070VXBFSfz+EdN7x5bH9pLx5ZkGzks2v9y9tXp5a+sOXb5bt\nX9519/hTSK+390NIjzn71zv2KOf5jqw7v/10d7q3gitC6u+nD+1un423n0Osby/cNv/yFO3d\n9a+8RHJ9nkHzFdj97yG93N5PIV3nXJa+c8euf+U6cXX97fH61y+//Xh3OreCK0Lq76eQLp9A\nLK9f+F/d21lXT48aH5rn2u3eImn+uFofjrcnLHy49ZvO7f0U0qnenD8WW+1ub/K4Y7e/0n3S\n3+V+fd8/qPtwd7qHhwtCmlXNZxpOENIsUvucndNh2X0CHYpFSLN4PFTw9vA0ikRIs7j/WwQe\n+nKCkOZRfzWPk1U8Y80LQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAAB\nQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUIC\nBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQI\nCRAgJECgd0j1OqXl7noj5Ijg+iZQV6mxutwIISG4vgls0vZc07ZatjdCSAiubwLV5Q2P1eJI\nSEDfBG7t1MslIQF9E1ik+vbSkpAQXt8Etml9femYloSE6HonsLnXs0uEhOj6J3BY3V46rgkJ\nwZEAIDBSSOnZOCOQJjT3sWZvgjPERRhJ+m8yXEMLIZWLkDJCSOUipIz0f2bDP38IzUUYCSFl\npP8XZAlpboSUkd5n6HB54veYI/A7QsrIgC/Ips3YI/ArQsrIgDO0TYexR+A3hJQRHrUrFyFl\nhJDKRUgZIaRyEVJGCKlchJQRQioXIWWEkMpFSBkhpHIRUkYIqVyElBFCKhchZYSQykVIGSGk\nchFSRgipXISUEUIqFyFlhJDKRUgZIaRyEVJGCKlchJQRQioXIWWEkMpFSBkhpHIRUkYIqVyE\nlBFCKhchZYSQykVIGSGkchFSRgipXISUEUIqFyFlhJDKRUgZIaRyEVJGCKlchJQRQioXIWWE\nkMpFSBkhpHIRUkYIqVyElBFCKhchZYSQykVIGSGkchFSRgipXISUEUIqFyFlhJDKRUgZIaRy\nEVJGCKlchJQRQioXIWWEkMpFSBkhpHIRUkYIqVyElBFCKhchZYSQykVIGSGkchFSRgipXISU\nEUIqFyFlhJDKRUgZIaRyEVJGCKlchJQRQioXIWWEkMpFSBkhpHIRUkYIqVyElBFCKhchZYSQ\nykVIGSGkchFSRgipXISUEUIqFyFlhJDKRUgZIaRyEVJGCKlchJQRQioXIWWEkMpFSBkhpHIR\nUkYIqVyElBFCKhchZYSQykVIGSGkchFSRgipXISUEUIqFyFlhJDKRUgZIaRyEVJGCKlchJSR\n/mdo/7VKjdVmP9YI/IqQMtL3DNWL9LAcZQQMhJSRvmdok6rvQ/vScVelzRgjYCCkjPQ9Q1U6\n3F8+pGqMETAQUkb6nqGUfvqNbAQMhJQR3iOVi5AyMuBzpN2xfYnPkeZCSBnpfYaWT4/aLepR\nRuB3hJSRAV9H2rRfR6pWX3wdaR6ElBGe2VAuQsrISGcoPRtnRK7SdAgpH7xHUptuvQkpI4Sk\nRkghEZIaIYXU/5kN//xpULCLQEgh9T1DW0L6ASGF1PsMHarf//GEYESZCCmk/mfo8PsTgxQj\nikRIIQ04Q9un562ONKJEhBQSj9qpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhoh\nhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJI\nhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURI\naoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRG\nSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQU\nEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCER\nkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGp\nEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhohhURIaoQUEiGpEVJIhKRGSCERkhoh\nhURIaoQUEiGpEVJIhKRGSCERkhohhdT7DB3Xqfo6nbaLVG1GGlEmQgqp7xmqq3S2/Wr+m5aj\njCgUIYXU9wxt0vn90KZK6/pUty/rRxSKkELqe4aq9g1TqttfqjFGFIqQQup7hlJ6/Pf2S+fV\nT/rfuxIRUkhD3yM1/615j/SEkEIa+jnSpr6+rB9RKEIKiUft1AgpJL6OpEZIIfHMBjVCComQ\n1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmN\nkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1Agp\nJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIi\nJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJSI6SQCEmNkEIiJDVCComQ1AgpJEJS\nI6SQns/Q4us49ogACCmk5zOUUhqjpWAXgZBCej5D9fd6jJaCXQRCCun1DO2/FuqWgl0EQgrp\nwxk6VOf3S9tRR3hGSCG9n6HdMjWWI45wjZBCejlD9df53dFiV59rWo00wjtCCqlzhvbNgw2b\nw+UVsnMX7CIQUkidryOd3xlt69srqjFGBEBIIXW+jrTajT0iAEIKqfN1pPFHBEBIIXXOUL1p\nPp6rNtqigl0EQgrp+Qwdq/YRhpQq6XMbgl0EQgrp+Qwt07p5X1RvdA99v44IgJBC6j5p9fUF\n+YgACCmk5zNUpcsnRzUhDUBIIT2foU1a7s+/7JdpM9aIAAgppM4ZujzLTvk8u7cR/hFSSN0z\n9L1qMhI+8/t9hHuEFBLfs0GNkEIiJDVCComQ1AgppM4Zav6Z+cVoI/wjpJCez9BXSoQ0GCGF\n1P2CrPjxuvcRARBSSB+fIjTeiAAIKaTnM7RKo/yLpGAXgZBC6v4zivYpQmOOCICQQnr5lsU8\n2DAYIYVESGqEFBJfkFUjpJAISY2QQuqeod2q+ahupf1xFMEuAiGF9P7vkZrvDck3P+mPkEJ6\nPkPbtGz/lfk2rccaEQAhhfT6PRuu35BrrBEBEFJIr08RIqShCCmk7jfRv7xHOqTFWCMCIKSQ\nPnyOtBM/CzzYRSCkkDpnaMV3ERqOkEJ6/zpSWn2POcI9QgqJZzaoEVJIhKRGSCH1PkOXn6XU\nfLuUpfGhYLCLQEgh9f1nFO3PUqqrf3lwIthFIKSQ+oa0Tqv6/J/18dzU+vdvuh/sIhBSSB/O\n0H75Dz9nLDXf3yFdvslD/ftPQA92EQgppE9nqP6HJ62277Sq9PSbl1eP9M9te0oTmm67J5tE\nSKaPZ+ifPrQ7NN9R8tC8XP/+SVIOF2HKnXM4iZBMn87Q9vcP1VqHVG0Op1V1Lmm3SLu/jpga\nIQ0cNfcFzN7nBxu+7DfcVf/613O4CIQ0cNTcFzB7n0Ja/NtzVr/X7ffcX30Z/542h4tASANH\nzX0BsxfkmQ2ENHDU3Bcwe4Qk3zmHkwjJ9MMXZBtjjJgLIQ0cNfcFzB4hyXfO4SRCMnV/Yl/V\nPI69r/z9wz5CGjhq7guYve5P7Gu/vno6pH94jlC/EXMhpIGj5r6A2fv4g8b8fRchQho4au4L\nmL3u97W7vUdy912ECGngqLkvYPaez9AmtZ8jefwuQoQ0cNTcFzB779/7++zXf140bMRMCGng\nqLkvYPa6Z+i7/S5Cvz4FdeiIeRDSwFFzX8Ds8cwG+c45nERIJkKS75zDSYRkCvKDxghp4Ki5\nL2D2gvygMUIaOGruC5i9ID9ojJAGjpr7AmYvyA8aI6SBo+a+gNkL8oPGCGngqLkvYPaC/KAx\nQho4au4LmL0gP2iMkAaOmvsCZi/IDxojpIGj5r6A2Qvyg8YIaeCouS9g9nhmg3znHE4iJNPz\nGVppn/X9acRcCGngqLkvYPY+/gvZ8UbMhZAGjpr7Ambv9eHvkUfMhZAGjpr7Ambv+QzVq+V+\n5BFzIaSBo+a+gNnr+xP7eo6YCyENHDX3BcweIcl3zuEkQjLx8Ld85xxOIiQTIcl3zuEkQjKN\n8z0hP46YEyENHDX3BcxeN6RRcsrhIhDSwFFzX8DsEZJ85xxOIiQTIcl3zuEkQjIRknznHE4i\nJBMhyXfO4SRCMhGSfOccTiIk0yOkUX7s5fOIORHSwFFzX8DsEZJ85xxOIiQTz2yQ75zDSYRk\nIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5\nzjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmc\nREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgm\nQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLv\nnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJ\nhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQi\nJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnO\nOZxESCZCku+cw0mEZBp+hpJ1EzlcBEIaOGruC5g9QpLvnMNJhGTqe4ZS1xgjlAhp4Ki5L2D2\n+p6hfUVIP+ycw0mEZOp9hupVWh7bW/h0E/9c2UQIaeCouS9g9gacoe+Uvk98jvS2cw4nEZJp\nyBk6LtOqJqTXnXM4iZBMw87QV6p2hPSycw4nEZJp4Bk6LOzPgXK4CIQ0cNTcFzB7g8/QmpBe\nds7hJEIy8RQh+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+c\nw0mEZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mE\nZCIk+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk\n+c45nERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45\nnERIJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERI\nJkKS75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nERIJkKS\n75zDSYRkIiT5zjmcREgmQpLvnMNJhGQiJPnOOZxESCZCku+cw0mEZCIk+c45nPRfmtDcu9IP\nIcl3zuEk3vmZCEm+CA4nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDk\ni+BwEiGZCEm+CA4nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+Bw\nEiGZCEm+CA4nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+BwEiGZ\nCEm+CA4nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+BwEiGZCEm+\nCA4nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+BwEiGZCEm+CA4n\nEZKJkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+BwEiGZCEm+CA4nEZKJ\nkOSL4HASIZkISb4IDicRkomQ5IvgcBIhmQhJvggOJxGSiZDki+BwEiGZet/tep3Scne9kV9v\nJYcz43LnXB5UtJDqKjVWlxshpKdFcDiJkEx97/Ymbc81batleyOE9LQIDicRkqnv3a4ub3is\nFkdC6i6Cw0mEZOp7t2/t1Mvlp5DSswF3T8Xlzrk8qGghLVJ9e2nJe6TOIjicREimvnd7m9bX\nl45pSUjPi+BwEiGZet/tzb2enfHRWw5nxuXOuTyocCGdDqvbS8c1IT0tgsNJhGTimQ3yRXA4\niZBMhCRfBIeTCMlESPJFcDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBM\nhCRfBIeTCMlESPJFcDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBMhCRf\nBIeTCMlESPJFcDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBMhCRfBIeT\nCMlESPJFcDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBMhCRfBIeTCMlE\nSPJFcDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBMhCRfBIeTCMlESPJF\ncDiJkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJREjyRXA4iZBMhCRfBIeTCMlESPJFcDiJ\nkEyEJF8Eh5MIyURI8kVwOImQTIQkXwSHkwjJNGdIaUITLoLDSZOOms4kWz7BCKeL4HCSz1GE\n1OOceRzl8qAI6e8jXF4dQipmFCH1OGceR7k8KEL6+wiXV4eQihlFSD3OmcdRLg+KkP4+wuXV\nIaRiRhFSj3PmcZTLgyKkv49weXUIqZhRhNTjnHkc5fKgCOnvI1xeHUIqZhQh9ThnHke5PChC\n+vsIl1eHkIoZRUg9zpnHUS4PipD+PsLl1SGkYkYRUo9z5nGUy4MipL+PcHl1CKmYUYTU45x5\nHOXyoAjp7yNcXh1CKmYUIfU4Zx5HuTwoQvr7CJdXh5CKGUVIPc6Zx1EuD4qQ/j7C5dUhpGJG\nEVKPc+ZxlMuDIqS/j3B5dQipmFGE1OOceRzl8qAI6e8jXF4dQipmFCH1OGceR7k8KEL6+wiX\nV4eQihlFSD3OmcdRLg+KkP4+wuXVIaRiRhFSj3PmcZTLgyKkv49weXUIqZhRhNTjnHkc5fKg\nCOnvI1xeHUIqZhQh9ThnHke5PChC+vsIl1eHkIoZRUg9zpnHUS4PipD+PsLl1SGkYkYRUo9z\n5nGUy4MipL+PcHl1CKmYUYTU45x5HOXyoAjp7yNcXh1CKmZUJiHtv1apsdrse45weXUIqZhR\nWYRUL9LDst8Il1eHkIoZlUVIm1R9H9qXjrsqbXqNcHl1CKmYUVmEVKXD/eVDqnqNcHl1CKmY\nUVmElNJPv7n+yZOfbwOYUc/d/7zMPd/uD++RAP8GfI60O7YvmZ8jAf71fve2fHoXuaiVdwko\nz4CvI23aryNVqy/j60iAfxM8swHwj5AAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIE\nCAkQICRAgJAAAUICBAgJECAkQICQAIE5Q5rpmzABF9JlVt5YvrNdjnJ5UKWeP0Iqd5TLgyr1\n/BFSuaNcHlSp54+Qyh3l8qBKPX+EVO4olwdV6vkjpHJHuTyoUs8fIZU7yuVBlXr+CKncUS4P\nqtTzR0jljnJ5UKWeP0Iqd5TLgyr1/BFSuaNcHlSp54+Qyh3l8qBKPX88+xsQICRAgJAAAUIC\nBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBGYLaXudPML3M/9p1GlT\npWpTjznqNMkBNSY5ltZEBzTldRph+eYK6XC9+4fxL9Nt1GnZTlqMOOo0yQE1JjmW1kQHNOV1\nGmP5ZgrpUN2PZTXVqH2qDs3v9uOOG/2AGtMcS2uaA5ryOo2yfPOEtE3L67Fs09dUozZpd/7v\n98gDRz+g1jTH0prmgCa8TuMs3zwhpc3pfizbqUat0vE0/v9gRz+g1jTH0prmgCa8TuMs3zwh\nHU6Ps7Zbnz+znGJU95exjH5ArWmOpTXNAU14ncZZvtketbsfS2s5waipQhr/gCYOaYoDakx3\nnUZYvrlDSun7dKo3o374MGlIUxzQpCFNc0CXSR9+GXeS7tDmDumiHvWB3ElDuhj3gKY9ltbY\nB9SYPKQLzaFNGtLzo/bd8yQ/ax9GVaNeoJevSIy84uMeyycTjJrkOn24bcmkQCFdHg06jvRo\n0LQhjXssn0wX0gTHVnpIncG3//00TwUZdyOuo77ar0/s0rgPQE1xQFMdS2uaA2pMd51GWL65\nQ9o056u+fA1u5FHTPBtgigOa9JkN0xxQY7rrNMLyzR1SXbUfEo36f9bbu+7FFI/jTnFAp4mO\npTXRAZ2mvE4jLN/cIZ3/h1ClxSSPFbejxv/K4gQHdJrqWO6jxj+g05TXaYTl498jAQKEBAgQ\nEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIg\nQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoRUmn/50ekT\n/AhydHHGS0NIWeKMF4yQ8sEZLxgh5YMznqOU6kVanV/YLlJ1+anbu2VKy93pFsmmSpv2xctv\n2//uVun608Avv729BSZASDlK6dzEOYnzf8+W5z/Zti+l7TWSZfObVTekr8tf2Vx/+3gLTICQ\ncnSOpz7/smt+qZfp/G6lSofT6TstLpF8p+pwOlTdkFL6bl6Trr99vAUmQEg5Smnf/LJKTU51\n80FeSrv765pXNK/fvXxod3/1tSs+rJsQIeXo2kW6OX9OdP5I7nC4ve7++m5Ix93X8h7S4y0w\nAULK0XtIp6/zB3KpOv4W0vL+l7tvgQkQUo6eQnnYbRa3z5E+hrROi+3u+Ajp/haYACHl6BrC\n6vXTnFs5l1fs7yE9XuqEdHptEaPhPOfouv7tg3OnbfNgw+LykNz1PdLu8ajdIm2bB/bS5RGK\nw+NzpMdbYAKElKPb+5HLZz3N5znfl89/9tfXtV9gWt+/XLS6Prrw/Hceb4EJEFKO7h+QbRfn\nXtrHC9rnKezvr/u6PbOheWl9eWnd/I3d5cHy57fABAgJECAkQICQAAFCAgQICRAgJECAkAAB\nQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUIC\nBAgJECAkQICQAAFCAgQICRAgJECAkACB/wF0l47XkefFWQAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"Histogram of residuals\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDT09PZ2dnh4eHp6enw8PD///8uNL8wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2d60LizBZEw0VUjorv/7SHSyAEid+AC63srPVDBaW6enfV\nIMho8ykiP6b5awMiFbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBY\nJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTA\nIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCea\npul/1F1xztOvmHmeNc1NK12Yve79u1v8gH8ZFLdaJrV3dxP/VKS32a9M7Hm79riLdDkoizQZ\n/qlIv5SHedO833aLtCJdXmeRJsOXIn3/Rb/k5XECj92KRZosQ/dIm+fF9qPl6+Ga5vhl66fd\nd1/r9iYf20uLl7Nbfsyb1faj1+X24/nq46j3Mm/mb5+fL7Nm8dZfvqf39f6l1dusZs3sIHdu\nrLvFx9Osmb9c3qke31/auVRp6V93vubVPXw7qL6hstTe3U0MFOlj1qZi0cvHov14ub/FW/sl\n3S3n+xscv6p5O1x7uPyxOl134lzvvK+fPb2jmbcLYyfPX42c7+jSzheVA/3rzte8vofvBnVh\nqCy1d3cTA0Xa3lFs/43dbJPwcpaP5TGShybNThePt2x2N3vZJmjz+bnq9aNpZme3bOnpXSvS\nXu+0zuzC2MnzVyNnO/pi54vKgf5152te38N3g7owVJbau7uJ5pz2isPb3Xc1m+2dwum6z/X2\n/ctm+83M9v32u7HXbbB272bdLXeJ3T1n8NFT2l67TXMzf9+/69a+0Lvyrd1e79CEzSGy14y1\nRtaz5lqRvtj5sr3Tct11vTWv7+GbQV0aKkvt3d3EQJF2ETg9FDrm4+n4D/hq/yz18hD/XWyO\nt1xfSB/evvXedV9woXelSHu97Tqbw+XldWNHI+urRfpi58v2DvSu6615fQ/fDOrSUFlq7+4m\nBor0fLiijUj3qX26Pj/2V8yOKbn89PYLXleL5jy5V59cv9C7UqRN3+LsurHmq5Hzi1fs9FUO\n9K7rrXl9D98PqmeoLLV3dxPdUfdDsDrm6OPLp44fDeX383V+1szvi9T76GuRju/P9K4Y+48i\nXbPTU2k5v66/5jdF+n5QFmkyDBXpc/N6eKJq0fvU6R5kdvUeaX9x963e/Onl/V+KdK43VKRZ\n/+qvxr4v0lU7PZUryv01vyvSwKAup1uT2ru7icEi7dj/lKe7bvmfj5H2n5231/9nkZb/+Rjp\n+GUXD2cujF0xsivoW3N8suGrnZ7KFeX+mt8W6fqgeobKUnt3NzFQpPnp8cnxrmIz+Kxdc/07\nqv++R/rvZ+3273frvO3fLa4Z2z/F1nv6cPt+tX/d27mBnp2+yoHedb01vyvS1UFdGipL7d3d\nxECRthlffOwfSu9eWbALxO796Uebh5eWDvz4ZrH/4nU/x9eKdKk3UKRunbdLY11zzo089S5e\ns9NTaW/1dcvHNb8r0vVB+XOkyTFQpNNj6P1jiKfjB8fkH16ivW6/5KIk7Y/1dz+9fPu+SBd6\nQ0Vq19lHtG+s/ZJ2xWV78aP98sPFq3bOVY7r9JR7a35TpOuDujBUltq7u4mhIh2+7V+0P/hf\nnrrzNDt71vh991q79ZeS7K6ePb1/HF+wcEX9SE9vqEifm9X2O6jl+niTztjxS/ov+jv4ej19\nU3fFzrlKV9jzLZ+v+V2Rrg/qwlBVau/ut9mcP9KQKWGRCJrD62feF/0X0Ml0sEgE3VMFX56e\nlmlgkQhO/4Pg8IhcpodFQtg8756Wml2++lMmg0USAbBIIgAWSQTAIokAWCQRAIskAmCRRAAs\nkgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJg\nkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEA\niyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokA\nWCQRgLuL9PHUzJ4/P1/mzWxFGhIZI/cWaTNrtrw87942C9SSyPi4t0irZns/tJo1T5vPzf5j\nkSlzb5Fm+xs2zWb/bob5ERkl9xapabq3x3e9T59xv7sx0kBUdFOYn94j7d5uvr9HmtghNP9D\noIoU5aYwP32MtNq0H/NLjJSs6Ga5KcwvPGs3sUPIim6Wm8L8ws+RJnYIWdHNclOYX5jQxA4h\nK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZ\nbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpj\nkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiy\nopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvl\npjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAW\niSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYr\nulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrullu\nCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmOR\naLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKi\nm+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+Wm\nMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJ\nJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6\nWW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4K\nY5FosqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5Fo\nsqKb5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb\n5aYwFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYw\nFokmK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokm\nK7pZbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZ\nbgpjkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpj\nkWiyopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiy\nopvlpjAWiSYrulluCmORaLKim+WmMBaJJiu6WW4KY5FosqKb5aYwFokmK7pZbgpjkWiyopvl\npjB3T2izmm3fPs+bZvH6oCXGSVZ0s9wU5t4Jfcya5nOzfbNj8ZAlRkpWdLPcFObeCT01y832\nzdPHtlNPzeoRS4yUrOhmuSnMvRNqmk37ZvtdXjN7xBIjJSu6WW4Kc3+Rtm9mzdmFi0+fcb+7\nf/XCALmJim6Wm8Lc/63d++fn8+7N7h7p2wdJjy9SVFh0M0nundB7M1u9fy5n2yat5836EUv8\nM1lh0c0kuXtC61n3PdHzY5b4V7LCoptJ8oMJvT7Ndy1aPn88bIl/IyssupkkFV7ZkBUW3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwSi0SHRTeTxCLRYdHNJLFIdFh0M0ksEh0W3UwS\ni0SHRTeTxCLRYdHNJLFIdFgoNxCQG2Y2hbFIdFgoN1kyzGwKY5HosIQ1AJJhZlMYi0SHJawB\nkAwzm8JYJDosYQ2AZJjZFMYi0WEJawAkw8ymMBaJDktYAyAZZjaFsUh0WMIaAMkwsymMRaLD\nEtYASIaZTWEsEh2WsAZAMsxsCmOR6LCENQCSYWZTGItEhyWsAZAMM5vCWCQ6LGENgGSY2RTG\nItFhCWsAJMPMpjAWiQ5LWAMgGWY2hbFIdFjCGgDJMLMpjEWiwxLWAEiGmU1hLBIdlrAGQDLM\nbApjkeiwhDUAkmFmUxiLRIclrAGQDDObwlgkOixhDYBkmNkUxiLRYQlrACTDzKYwFokOS1gD\nIBlmNoWxSHRYwhoAyTCzKcz5hObPH49e4iFkhSWsAZAMM5vCnE+oaZpHdMki3eUmS4aZTWHO\nJ7R5fXpElyzSXW6yZJjZFOZyQm/Pc7pLFukuN1kyzGwKc2VC77Pt/dLLQ5dgyQpLWAMgGWY2\nhfk6ofVi/8vXFw9cAiYrLGENgGSY2RTmYkKb5+3d0Xy92bZp+aAleLLCEtYASIaZTWF6E3rb\nPdmwej98ApudRbrLTZYMM5vC9H6OtL0zetkcPzF7xBIPISssYQ2AZJjZFKb3c6Tl+tFLPISs\nsIQ1AJJhZlOY3s+RHr/EYxaICktYAyAZZjaF6U1os9p9PzdbsY2ySHe5yZJhZlOY8wl9zPbP\nMDTNDH1tg0W6y02WDDObwpxPaNE87e6LNivuqe/LJR5CVljCGgDJMLMpTP9Fq5cf4Es8hKyw\nhDUAkmFmU5jzCc2aw4OjjUUKcJMlw8ymMOcTWjWLt+27t0WzetQSDyErLGENgGSY2RSmN6HD\nq+zI19l9WeIRZIUlrAGQDDObwvQn9Lrc1Qh85ffXJR5AVljCGgDJMLMpjL+zgQ5LWAMgGWY2\nhbFIdFjCGgDJMLMpjEWiwxLWAEiGmU1hehPa/TfzAw9b4hFkhSWsAZAMM5vCnE/ouWksUoyb\nLBlmNoXp/0AWfr7u6xIPISssYQ2AZJjZFObqS4Qet8RjFogKS1gDIBlmNoU5n9Cyecj/SLJI\nd7nJkmFmU5j+f6PYv0TokUs8hKywhDUAkmFmU5iLX1nskw0xbrJkmNkUxiLRYQlrACTDzKYw\n/kCWDktYAyAZZjaFsUh0WMIaAMkwsylMf0Lr5e67uiX75ygs0l1usmSY2RTm6/9H2v1uSH/5\nyd+7yZJhZlOY8wm9NIv9/zJ/aZ4etcRDyApLWAMgGWY2hbn8nQ3tL+R61BIPISssYQ2AZJjZ\nFObyJUIWKcVNlgwzm8KcT2je3iO9N/NHLfEQssIS1gBIhplNYa48RlrDrwK3SHe5yZJhZlOY\n3oSW/hahHDdZMsxsCvP150jN8vWRSzyArLCENQCSYWZTGF/ZQIclrAGQDDObwlgkOixhDYBk\nmNkU5ucT+s/nyi3SXW6yZJjZFObn/43CIj3GTZYMM5vC3Fukps+/LvEQssIS1gBIhplNYa5M\n6G3xD39n7G1mkR7rJkuGmU1hrk1o8y8vWt0sm8X+ReJXW/TPLSPICktYAyAZZjaFuTqhf8v+\na9O8/ssXW6S73GTJMLMpzLUJvTSzf7rtx6JZbizSg9xkyTCzKcz1Jxue//HWz81sbZEe4yZL\nhplNYa4Vaf7vr1l9n//3YyCLdJebLBlmNoX58YSeLNJj3GTJMLMpjC8RosMS1gBIhplNYQZ+\nIEs+a22R7nKTJcPMpjAWiQ5LWAMgGWY2helN6Hm2/ty9aMH/2BfgJkuGmU1hzif03Lzv3783\n//AaofuWeAhZYQlrACTDzKYwl79FqP8BvsRDyApLWAMgGWY2hen/XrvjPZK/Rejv3WTJMLMp\nzPmEVs3+MZK/RSjCTZYMM5vC9Ca0aJ+vWz1uiUeQFZawBkAyzGwK05/Q6/63CK0fucQDyApL\nWAMgGWY2hfGVDXRYwhoAyTCzKYxFosMS1gBIhplNYb7+gkj/0FiGmywZZjaF+fpkw6d/aCzC\nTZYMM5vCnE/IPzSW5CZLhplNYfo/kPUPjeW4yZJhZlOYy5cIWaQUN1kyzGwKcz4h/9BYkpss\nGWY2hbnyGMmXCEW4yZJhZlOY3oT8Q2NBbrJkmNkU5uvPkfxDYxlusmSY2RTGVzbQYQlrACTD\nzKYw5xNasq/6vrbEYxaICktYAyAZZjaFufo/ZB+3xGMWiApLWAMgGWY2hbl8+vvBSzyErLCE\nNQCSYWZTmPMJbZaLtwcv8RCywhLWAEiGmU1hBn6v3aOWeAhZYQlrACTDzKYwFokOS1gDIBlm\nNoXx6W86LGENgGSY2RTGItFhCWsAJMPMpjCP+Z2QV5d43AJRYQlrACTDzKYw/SI9pE4W6S43\nWTLMbApjkeiwhDUAkmFmUxiLRIclrAGQDDObwlgkOixhDYBkmNkUxiLRYQlrACTDzKYwFokO\nS1gDIBlmNoXpivSQP3t5vsTDyApLWAMgGWY2hbFIdFjCGgDJMLMpjK9soMMS1gBIhplNYSwS\nHZawBkAyzGwKY5HosIQ1AJJhZlMYi0SHJawBkAwzm8JYJDosYQ2AZJjZFMYi0WEJawAkw8ym\nMBaJDktYAyAZZjaFsUh0WMIaAMkwsymMRaLDEtYASIaZTWEsEh2WsAZAMsxsCmOR6LCENQCS\nYWZTGItEhyWsAZAMM5vCWCQ6LGENgGSY2RTGItFhCWsAJMPMpjAWiQ5LWAMgGWY2hbFIdFjC\nGgDJMLMpjEWiwxLWAEiGmU1hLBIdlrAGQDLMbApjkeiwhDUAkmFmUxiLRIclrAGQDDObwlgk\nOixhDYBkmNkUxiLRYQlrACTDzKYwFokOS1gDIBlmNoWxSHRYwhoAyTCzKYxFosMS1gBIhplN\nYSwSHZawBkAyzGwKY5HosIQ1AJJhZlMYi0SHJawBkAwzm8JYJDosYQ2AZJjZFMYi0WEJawAk\nw8ymMBaJDktYAyAZZjaFsUh0WMIaAMkwsymMRaLDEtYASIaZTWEsEh2WsAZAMsxsCmOR6LCE\nNQCSYWZTGItEhyWsAZAMM5vCWCQ6LGENgGSY2RTGItFhCWsAJMPMpjAWiQ5LWAMgGWY2hbFI\ndFjCGgDJMLMpjEWiwxLWAEiGmU1hLBIdlrAGQDLMbApjkeiwhDUAkmFmUxiLRIclrAGQDDOb\nwlgkOixhDYBkmNkUxiLRYQlrACTDzKYwFokOS1gDIBlmNoWxSHRYwhoAyTCzKYxFosMS1gBI\nhplNYSwSHZawBkAyzGwKY5HosIQ1AJJhZlMYi0SHJawBkAwzm8JYJDosYQ2AZJjZFMYi0WEJ\nawAkw8ymMBaJDktYAyAZZjaFsUh0WMIaAMkwsymMRaLDEtYASIaZTWEsEh2WsAZAMsxsCmOR\n6LCENQCSYWZTGItEhyWsAZAMM5vCWCQ6LGENgGSY2RTGItFhCWsAJMPMpjAWiQ5LWAMgGWY2\nhbFIdFjCGgDJMLMpjEWiwxLWAEiGmU1hLBIdlrAGQDLMbApjkeiwhDUAkmFmUxiLRIclrAGQ\nDDObwlgkOixhDYBkmNkUxiLRYQlrACTDzKYwFokOS1gDIBlmNoWxSHRYwhoAyTCzKYxFosMS\n1gBIhplNYSwSHZawBkAyzGwKY5HosIQ1AJJhZlMYi0SHJawBkAwzm8JYJDosYQ2AZJjZFMYi\n0WEJawAkw8ymMBaJDktYAyAZZjaFsUh0WMIaAMkwsymMRaLDEtYASIaZTWEsEh2WsAZAMsxs\nCmOR6LCENQCSYWZTGItEhyWsAZAMM5vCWCQ6LGENgGSY2RTGItFhCWsAJMPMpjAWiQ5LWAMg\nGWY2hbFIdFjCGgDJMLMpjEWiwxLWAEiGmU1hLBIdlrAGQDLMbApjkeiwhDUAkmFmUxiLRIcl\nrAGQDDObwlgkOixhDYBkmNkUxiLRYQlrACTDzKYwFokOS1gDIBlmNoW5e0Kbp6ZZrFuRb1Us\n0l1usmSY2RTm3gltZs2O5UHEIvFusmSY2RTm3gmtmpdtm15mi72IReLdZMkwsynMvROaHW74\nMZt/WKSHuMmSYWZTmHsndOzOZrG4VqTmnB/Y+zcvUFggIDdZMtBJQTBuUO71NG82x48WVe6R\nlBmWyTopxg3KvZ5emqf2o49mYZHKy2SdFOMG5W5Pq1N71v9xX2uRCshknRTjBuV+T+/L40cf\nTxapukzWSTFuUHxlQ3c8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PM\nsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJ\nMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAs\nUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzK\nDMtknRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2Sd\nFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbF\nInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcej\nzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzW\nSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5Q\nLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAsUnc8\nygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtk\nnRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMG\nxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblAsUnc8ygzLZJ0U4wbFInXH\no8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFuUCxSdzzKDMtknRTjBsUidcejzLBM\n1kkxblAsUnc8ygzLZJ0U4wbFInXHo8ywTNZJMW5QLFJ3PMoMy2SdFOMGxSJ1x6PMsEzWSTFu\nUCxSdzzKDMtknRTjBsUidcejzLBM1kkxblD+skgNBXQ8ygzLQFGA3EAwm2q3RorduAQ01rTM\nlZSBogC5gWSYTbVbI8VuXMIijUcGigLkBpJhNtVujRS7cQmLNB4ZKAqQG0iG2VS7NVLsxiUs\n0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFkmE21WyPFblzCIo1HBooC5AaSYTbVbo0Uu3EJ\nizQeGSgKkBtIhtlUuzVS7MYlLNJ4ZKAoQG4gGWZT7dZIsRuXsEjjkYGiALmBZJhNtVsjxW5c\nwiKNRwaKAuQGkmE21W6NFLtxCYs0HhkoCpAbSIbZVLs1UuzGJSzSeGSgKEBuIBlmU+3WSLEb\nl7BI45GBogC5gWSYTbVbI8VuXMIijUcGigLkBpJhNtVujRS7cQmLNB4ZKAqQG0iG2VS7NVLs\nxiUs0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFkmE21WyPFblzCIo1HBooC5AaSYTbVbo0U\nu3EJizQeGSgKkBtIhtlUuzVS7MYlLNJ4ZKAoQG4gGWZT7dZIsRuXsEjjkYGiALmBZJhNtVsj\nxW5cwiKNRwaKAuQGkmE21W6NFLtxCYs0HhkoCpAbSIbZVLs1UuzGJSzSeGSgKEBuIBlmU+3W\nSLEbl7BI45GBogC5gWSYTbVbI8VuXMIijUcGigLkBpJhNtVujRS7cQmLNB4ZKAqQG0iG2VS7\nNVLsxiUs0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFkmE21WyPFblzCIo1HBooC5AaSYTbV\nbo0Uu3EJizQeGSgKkBtIhtlUuzVS7MYlLNJ4ZKAoQG4gGWZT7dZIsRuXsEjjkYGiALmBZJhN\ntVsjxW5cwiKNRwaKAuQGkmE21W6NFLtxCYs0HhkoCpAbSIbZVLs1UuzGJSzSeGSgKEBuIBlm\nU+3WSLEbl7BI45GBogC5gWSYTbVbI8VuXMIijUcGigLkBpJhNtVujRS7cQmLNB4ZKAqQG0iG\n2VmM/Z8AAAiPSURBVFS7NVLsxiUs0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFkmE21WyPF\nblzCIo1HBooC5AaSYTbVbo0Uu3EJizQeGSgKkBtIhtlUuzVS7MYlLNJ4ZKAoQG4gGWZT7dZI\nsRuXsEjjkYGiALmBZJhNtVsjxW5cwiKNRwaKAuQGkmE21W6NFLtxCYs0HhkoCpAbSIbZVLs1\nUuzGJSzSeGSgKEBuIBlmU+3WSLEbl7BI45GBogC5gWSYTbVbI8VuXMIijUcGigLkBpJhNtVu\njRS7cQmLNB4ZKAqQG0iG2VS7NVLsxiUs0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFkmE21\nWyPFblzCIo1HBooC5AaSYTbVbo0Uu3EJizQeGSgKkBtIhtlUuzVS7MYlLNJ4ZKAoQG4gGWZT\n7dZIsRuXsEjjkYGiALmBZJhNtVsjxW5cwiKNRwaKAuQGkmE21W6NFLtxCYs0HhkoCpAbSIbZ\nVLs1UuzGJSzSeGSgKEBuIBlmU+3WSLEbl7BI45GBogC5gWSYTbVbI8VuXMIijUcGigLkBpJh\nNtVujRS7cQmLNB4ZKAqQG0iG2VS7NVLsxiUs0nhkoChAbiAZZlPt1kixG5ewSOORgaIAuYFk\nmE21WyPFblzCIo1HBooC5AaSYTbVbo0Uu3EJizQeGSgKkBtIhtlUu7W7b/n2vGx2LFdvdy5h\nkcYjc3dMHnHiWZtqt3bn7TbzpmNx3xIWaTwyd8bkMSeetal2a3febtXMXt/3H32sZ83qriUs\n0nhk7ozJY048a1Pt1u683ax5P3383szuWsIijUfmzpg85sSzNtVu7d7bNUMX2mvOGNYQ+UPu\nzP71MN95uxvukUTq84PHSOuP/Uf/+RhJpD53370tzu4i5xvSksj4+MHPkVb7nyPNls//8XMk\nkfr8wisbROpjkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJ\nBMAiiQBYJBEAiyQCYJFEAP6ySH/0S5hk9PxhaIf40yL94dpXyLKjm2Gy3BywSCey7OhmmCw3\nByzSiSw7uhkmy80Bi3Qiy45uhslyc8Ainciyo5thstwcsEgnsuzoZpgsNwcs0oksO7oZJsvN\nAYt0IsuObobJcnPAIp3IsqObYbLcHLBIJ7Ls6GaYLDcHLNKJLDu6GSbLzQGLdCLLjm6GyXJz\nINGTyOiwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCR\nRAAskgjAnxbpZd7MVpu/dHDBS8i/K6tZ1mBS5rIjLjQH/nJCq/1fFpjlDOU95O8cLPaDmf+1\njRMpc9kRF5qWP5zQe/O02f1j9/R3Fvq8zzIC89bM3ndm3v7aSEvKXHbEhebIH05oeVg75pBe\nmkWGl1Wz3r59bZ7/2siBmLnsSAvNib83FDOTZhXiZdl8fO7+7V3+tZEDMXM5I87Q3xdp0yz+\n2kLLe8r5NFn/6sbMpSMnNCf+fEIv++9jQsgITFiRPqOs7IkKzYG/ntDHLOQbmD0ZgbFI/0FW\naA788YQ2s6j76IzAWKTvCQvNgT+Y0Pkfpl78/Q9Lzu1kBGZmkb4lIDRf+dMifcwXH7+//gV5\nRTo8a/eR8qzdZ8pcDkSE5it/OaF13HMvGYF53j+UXjervzZyImMue/JCc+APJ/SRN5KMwKS9\nsiFlLjsCQ3PgDyf01DTn31YlEOJlvh9LUGJC5vIZGZoDf2ioyZtJiJfN/tXff+3ijJC5fEaG\n5kCcIZExYpFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTA\nIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIA\nFkkEwCKJAFgkEQCLJAJgkYI4+zt0g3+S7u6/VXdxwzWhKSccYRC/WKR5M/gpuQdHmAlfpG9k\nLNLPcYSZWKSR4Qh/mabZzJvl9oOXeTN72V+1XjTNYv15DPRq1qz2Hx4u7t+ul037d84PF4+3\naDm7fJLdfuFzM3ve6jXN6vKGpz8MflpMfoYj/GWaZrkP9nKf5cX2mpfD37t/abO+2F1Y9ov0\nfPiSYx+6Wxw4u9zJNs3+VutFe8veDY9F6haTn+EIf5ltyjfbd+vdu82i2d47zJr3z8/XZn7I\n+msze/98n/WL1DSvu8807cXuFge6y2ey+4Ve2rezyxsetM8Wk5/hCH+ZpnnbvVs2uzptdt/k\nNc369LndJ3afX198a3f6dNur9YXm8XJP9m3/qY/TjS6X6i0mP8MR/jJtaJsj+8cwy/f34+dO\nn+8X6WP9vDgVqbvFgXOFk+z5rY+X+0v1FpOf4Qh/ma9F+nzefm/VzD6+K9Li9MX9W7ScKXxb\npP5SFgnEEf4yV7O7Xs2PD1yuFumpmb+sP7oinW5xTaG30EWRLr7QImE4wl+mDe3y4mHOKeuH\nT7ydot991CvS59f4d7c++/SXIvUunS0mP8MR/jJtaPfPl32+7J4VmB+ekmvvJtbdE2nz5mX3\nDFxzeOLgvXuM1N3iQHf5TPZakc6X2n1juPZZOwpH+MscQ3t41LN7uPJ6eFTz1n5u/5Ogp9NP\nfZbtkwTnX9Pd4nCTs8ud7LUidV84b3ZPiZ8tJj/DEf4yp9C+bLP8tH++YP9yg1Mrdk8IrE5P\nDTwdPnrafcX6dEdzukV7k+5yJ3v1W7vTF77N90U6W0x+hCMUAbBIIgAWSQTAIokAWCQRAIsk\nAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgk\nEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCRRAD+DwuyORqY78+YAAAAAElF\nTkSuQmCC",
      "text/plain": [
       "Plot with title \"Histogram of residuals.semistd\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDT09PZ2dnh4eHp6enw8PD///8uNL8wAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2d6UIi6xIEm0VUrgrv/7SXpaEBYQ6NAWQXET/cyS+ryBwE\nGW2WIvJnmmcbEKmARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJ\nAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJ\nBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEW6hqZp\njt/qPnDI20PMvI+aptdJJ2bPe//XJf7BsZHrL1eOlx28F1cV6Wv0kGW+r85OKdLpyBZJ/slV\nRXpQisZN893vEvcr0unXWST5J7+K9O8vepCX+wlYpN687OC9uHSLtHifrN6afm4/0uy+bP62\n/u5r3l7kZ/Xe5OPgkj/jZrZ663O6ens8+9npfYyb8ddy+TFqJl/Hxx/p/b59afUWs1Ez2sod\nGusu8fM2asYfpzequ9endk5VWs6PfCz9irzs4L24UKSfUZulyVGqJu3b080lvtov6S453lxg\n91XN1/aj2/d/ZvuP7TnUO+zr8khvZ+brxNje828jhxOd2vmlsuX8yCfSr8jLDt6LC0Va3VCs\n/mVerPLzcZCq6S6S2yaN9u/uLtmsL/axyt1iuZwd9aNpRocd3HKkd65IG739OaMTY3vPv40c\nTPTLzi+VLedHPpF+RV528F40h7Qf2L5cfye0WN0o7D+2nK9efyxW3wKtXq++G/tcRXv9atRd\ncp3Y9WMGP0dKq4+u0tyMvzevurNP9M58a7fR2zZhsQ36OWOtkfmoOVekX3Z+jbc/7vfIp9Kv\nyMsO3osLRVoHZ39XaJeqt90/4LPNo9TTbfzXYdtdcn4ivX35dfSq+4ITvTNF2uitzlls35+e\nN7YzMj9bpF92fo235ezIp9KvyMsO3osLRXrffqANVvepTaKXP5sPjHbZOv306gs+Z5PmMLln\nH1w/0TtTpMWxxdF5Y81vI4fvnrFzrLLl0sjLY63X42UH70UXkOPozHbZ/fn1qd1bl/K7/Bwf\nNPPfRTp663eRdq8P9M4Y+48inbNzpNLyz5EtkvyTS0VaLj63j3ZNjj61vwUZnb1F2ry7/lZv\n/Pbx3esWaXT6ye7d0fGHfxv7d5HO2jlSuaxskZYW6TouFmnN5qc83cem/3kfafPZcfvx/yzS\n9D/vI+2+7OTuzImxM0bWBf1qdg82/LZzpHJZ+Zf0K/Kyg/fiQpHG+/snu5uKxcVH7Zrz31H9\n9y3Sfz9qt3m9Pudr82pyztjmYb2jhw9Xr2ebZ8sdGjiyc6yy5ezIp9KvyMsO3osLRVplfPKz\nuQO+fmbBOkbr1/sfbW6fWnrhxzeTzRfPj3N8rkineheK1J3zdWqsa86hkbejd8/ZOVJpL3V+\nZH+OZJGu4tK3drt73pv7EG+7N3bJ3z5Fe95+yUlJ2icDrH8C+/XvIp3oXSpSe84m2MfG2i9p\nT5y27/60X75996ydQ5XdOWdHPpF+RV528F5cKtL2zsKk/cH/dN+dt9HBo8bf6+fazX+VZP3h\n0dv3z+4JC2fUdxzpXSrScjFbfd81ne8u0hnbfcnxk/62vj7339SdsXOo0hX23Mgn0i/Iyw7+\naBaH9zSkHBbpzjSb5+wsvyfHT6CTYlikO9M9VPDr4WkphEW6M/v/d7B9FECKYpHuzeJ9/WDW\n6PTZn1ILiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAW\nSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGw\nSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiA\nRRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkWrT\nQDx7jnjcUG2a/yEYk//CDdXGIj0IN1Qbi/Qg3FBtLNKDcEO1sUgPwg3VxiI9CDdUG4v0INxQ\nbSzSg3BDtbFID8IN1cYiPQg3VBuL9CDcUG0s0oNwQ7WxSA/CDdXGIj0IN1Qbi/Qg3FBtLNKD\ncEO1sUgPwg3VxiI9CDdUG4v0INxQbSzSg3BDtbFID8IN1cYiPQg3VBuL9CDcUG0s0oNwQ7Wx\nSA/CDdXGIj0IN1Qbi/Qg3FBtLNKDcEO1sUgPwg3VxiI9CDdUG4v0INxQbSzSg3BDtbFID8IN\n1cYiPQg3VBuL9CDcUG0s0oNwQ7WxSA/CDdXGIj0IN1Qbi/Qg3FBtLNKDcEO1sUgPwg3VxiI9\nCDdUG4v0INxQbSzSg3BDtbFID8IN1cYiPQg3VBuL9CDcUG0s0oNwQ7WxSA/CDdXGIj0IN1Qb\ni/Qg3FBtLNKDcEO1sUgPwg3VxiI9CDdUG4v0INxQbSzSg3BDtbFID8IN1cYiPQg3VBuL9CBu\n3tDPWzN6Xy4/xs1oRhoSFIv0IG7d0GLUrPh4X79sJqglAbFID+LWDc2a1e3QbNS8LZaLzdsS\niUV6ELduaLS5YNMsNq9GmB9hsUgP4tYNNU33cvfq6NMH3O5uiDQQkBuL9Bj+eou0frn49y3S\ni10JWdHNclOYv95Hmi3at/kjBkpWdLPcFOYBj9q92JWQFd0sN4V5wM+RXuxKyIpulpvCPGBD\nL3YlZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLR\nZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3\ny01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01h\nLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJN\nVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz\n3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTG\nItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFk\nRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfL\nTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEs\nEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1W\ndLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPc\nFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi\n0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRF\nN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tN\nYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwS\nTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0\ns9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wU\nxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLR\nZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFObmDS1mo9XL93HTTD7vdMQwyYpulpvC\n3Lqhn1HTLBerF2smdzlioGRFN8tNYW7d0FszXaxevP2sOvXWzO5xxEDJim6Wm8LcuqGmWbQv\nVt/lNaN7HDFQsqKb5aYwtxdp9WLUHLxz8ukDbnc3RLKim+WmMLd/a/e9XL6vX6xvkf55J+nF\nroSs6Ga5KcytG/puRrPv5XS0atJ83MzvccRAyYpulpvC3Lyh+aj73u39PkcMk6zoZrkpzB82\n9Pk2Xrdo+v5ztyOGSFZ0s9wUxmc20GRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTG\nItFkRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFk\nRTfLTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfL\nTWEsEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEs\nEk1WdLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1W\ndLPcFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPc\nFMYi0WRFN8tNYSwSTVZ0s9wUxiLRZEU3y01hLBJNVnSz3BTGItFkRTfLTWEsEk1WdLPcFMYi\n0WRFN8tNYQ43NH7/ufcRL0BWdLPcFOZwQ03T3KNLL3YlZEU3y01hDje0+Hy7R5de7ErIim6W\nm8KcbujrfUx36cWuhKzoZrkpzJkNfY9Wt0sfdz2iMlnRzXJTmN8bmk+aNZM7HlGarOhmuSnM\nyYYW76ubo/F8sWrT9E5HVCcrulluCnO0oa/1gw2z7+0nsN292JWQFd0sN4U5+jnS6sboY7H7\nxOgeR7wAWdHNclOYo58jTef3PuIFyIpulpvCHP0c6f5HvABZ0c1yU5ijDS1m6+/nRjO2US92\nJWRFN8tNYQ439DPaPMLQNCP0uQ0vdiVkRTfLTWEONzRp3ta3RYsZ99D36REvQFZ0s9wU5vhJ\nq6dv4Ee8AFnRzXJTmMMNjZrtnaOFRfoDWdHNclOYww3NmsnX6tXXpJnd64gXICu6WW4Kc7Sh\n7bPsyOfZ/TqiPlnRzXJTmOMNfU7XNQKf+f37iPJkRTfLTWH8nQ00WdHNclMYi0STFd0sN4Wx\nSDRZ0c1yU5ijDa3/m/mWux1Rn6zoZrkpzOGG3pvGIv2ZrOhmuSnM8Q9k4cfrfh/xAmRFN8tN\nYc4+Reh+R7wAWdHNclOYww1Nm7v8j6QXuxKyopvlpjDH/41i8xShex7xAmRFN8tNYU5+ZbEP\nNvyZrOhmuSmMRaLJim6Wm8L4A1marOhmuSmMRaLJim6Wm8Icb2g+XX9XN2X/HMWLXQlZ0c1y\nU5jf/x9p/bsh/eUnt5MV3Sw3hTnc0Ecz2fwv84/m7V5HvABZ0c1yU5jT39nQ/kKuex3xAmRF\nN8tNYU6fImSR/kpWdLPcFOZwQ+P2Fum7Gd/riBcgK7pZbgpz5j7SHH4W+ItdCVnRzXJTmKMN\nTf0tQn8nK7pZbgrz++dIzfTznkeUJyu6WW4K4zMbaLKim+WmMBaJJiu6WW4KY5FosqKb5aYw\n/jcKmqzoZrkpjEWiyYpulpvCnNnQ1wT9O2MW6ZnRzXJTmHMbWvik1T+QFd0sN4U5uyG/tfsD\nWdHNclOYcxv6aEb3PqIwWdHNclOY8w82vN/riBcgK7pZbgpzrkhj9jcXv9iVkBXdLDeF8Qey\nNFnRzXJTGItEkxXdLDeFufADWfKHsi92JWRFN8tNYSwSTVZ0s9wU5mhD76P56uXXyP/Y9wey\nopvlpjCHG3pvvjevvxv0OUIvdiVkRTfLTWFOf4vQ8Rv4ES9AVnSz3BTm+Pfa7W6R/C1Ct5MV\n3Sw3hTnc0KzZ3Efytwj9iazoZrkpzO/f/b1idr8j6pMV3Sw3hTne0OfmtwjN73lEebKim+Wm\nMD6zgSYrulluCmORaLKim+WmML9/QaR/aOxvZEU3y01hfj/YsPQPjf2JrOhmuSnMmV+i7x8a\n+xNZ0c1yU5jjH8j6h8b+TlZ0s9wU5vQpQhbpr2RFN8tNYQ435B8aI8iKbpabwpy5j+RThP5E\nVnSz3BTmaEP+oTGArOhmuSnM758j+YfG/kZWdLPcFMZnNtBkRTfLTWEONzS96Vnf//kQ34td\nCVnRzXJTmLP/Q7afgkU6Iiu6WW4Kc/rw99WXu/43Dr3YlZAV3Sw3hTnc0GI6+br2cl8ji3Se\nrOhmuSnMhd9r998XXEybyc/2Uudk7/EL8i7RQEBuoOhCQG6Y3RTm5iItl59N87lMuI+UFRbK\nTZYMs5vC/GVDP5NmurBId3KTJcPspjB/29B7M5pbpPu4yZJhdlOYP/5OyO/xf38faJFucpMl\nw+ymMMdFuqFObxbpPm6yZJjdFObPRbr+iLuRFZawBkAyzG4KY5HosIQ1AJJhdlMYi0SHJawB\nkAyzm8JYJDosYQ2AZJjdFMYi0WEJawAkw+ymMF2R7vBsmeMj7kZWWMIaAMkwuymMRaLDEtYA\nSIbZTWEq/A/ZrLCENQCSYXZTGItEhyWsAZAMs5vCWCQ6LGENgGSY3RTGItFhCWsAJMPspjAW\niQ5LWAMgGWY3hbFIdFjCGgDJMLspjEWiwxLWAEiG2U1hLBIdlrAGQDLMbgpjkeiwhDUAkmF2\nUxiLRIclrAGQDLObwlgkOixhDYBkmN0UxiLRYQlrACTD7KYwFokOS1gDIBlmN4WxSHRYwhoA\nyTC7KYxFosMS1gBIhtlNYSwSHZawBkAyzG4KY5HosIQ1AJJhdlMYi0SHJawBkAyzm8JYJDos\nYQ2AZJjdFMYi0WEJawAkw+ymMBaJDktYAyAZZjeFsUh0WMIaAMkwuymMRaLDEtYASIbZTWEs\nEh2WsAZAMsxuCmOR6LCENQCSYXZTGItEhyWsAZAMs5vCWCQ6LGENgGSY3RTGItFhCWsAJMPs\npjAWiQ5LWAMgGWY3hbFIdFjCGgDJMLspjEWiwxLWAEiG2U1hLBIdlrAGQDLMbgpjkeiwhDUA\nkmF2UxiLRIclrAGQDLObwlgkOixhDYBkmN0UxiLRYQlrACTD7KYwFokOS1gDIBlmN4WxSHRY\nwhoAyTC7KYxFosMS1gBIhtlNYSwSHZawBkAyzG4KY5HosIQ1AJJhdlMYi0SHJawBkAyzm8JY\nJDosYQ2AZJjdFMYi0WEJawAkw+ymMBaJDktYAyAZZjeFsUh0WMIaAMkwuymMRaLDEtYASIbZ\nTWEsEh2WsAZAMsxuCmOR6LCENQCSYXZTGItEhyWsAZAMs5vCWCQ6LGENgGSY3RTmmUVqKKLC\nEtYASIbZTWGeWiTmSg4LS1gDIBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDj\nQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6y\nZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjd\nFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYi\ndTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQ\nuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCb\nLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlm\nN4WxSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4Wx\nSJ0MNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0M\nNC7kJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7k\nJkuG2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG\n2U1hLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1h\nLFInA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFIn\nA40LucmSYXZTGIvUyUDjQm6yZJjdFMYidTLQuJCbLBlmN4WxSJ0MNC7kJkuG2U1hLFInA40L\nucmSYXZTmJs3tHhrmsm8FfmnikW6yU2WDLObwty6ocWoWTPdilgkfKowGWY3hbl1Q7PmY9Wm\nj9FkI2KR8KnCZJjdFObWDY22F/wZjX8s0j2mCpNhdlOYWze0685iMjlXpOaQixrMlYyFBQJy\nkyVzY0z+FYs/wLhBudXTuFns3ppUuUVS5rLMjTG5zzVeqUgfzVv71k8zsUjlZW6MyX2u8UpF\nWs727Zn/x22tRSogc2tM7nKNlyrS8nu6e+vnzSJVl7k5Jve4xmsV6e9HWKThyEBRiHKDYpGU\nuUYGikKUGxSLpMw1MlAUotygWCRlrpGBohDlBsUiKXONDBSFKDcoFkmZa2SgKES5QbFIylwj\nA0Uhyg2KRVLmGhkoClFuUCySMtfIQFGIcoNikZS5RgaKQpQbFIukzDUyUBSi3KBYJGWukYGi\nEOUGxSIpc40MFIUoNygWSZlrZKAoRLlBsUjKXCMDRSHKDYpFUuYaGSgKUW5QLJIy18hAUYhy\ng2KRlLlGBopClBsUi6TMNTJQFKLcoFgkZa6RgaIQ5QbFIilzjQwUhSg3KBZJmWtkoChEuUGx\nSMpcIwNFIcoNikVS5hoZKApRblAskjLXyEBRiHKDYpGUuUYGikKUGxSLpMw1MlAUotygWCRl\nrpGBohDlBsUiKXONDBSFKDcoFkmZa2SgKES5QbFIylwjA0Uhyg2KRVLmGhkoClFuUCySMtfI\nQFGIcoNikZS5RgaKQpQbFIukzDUyUBSi3KBYJGWukYGiEOUGxSIpc40MFIUoNygWSZlrZKAo\nRLlBsUjKXCMDRSHKDYpFUuYaGSgKUW5QLJIy18hAUYhyg2KRlLlGBopClBsUi6TMNTJQFKLc\noFgkZa6RgaIQ5QbFIilzjQwUhSg3KBZJmWtkoChEuUGxSMpcIwNFIcoNikVS5hoZKApRblAs\nkjLXyEBRiHKDYpGUuUYGikKUGxSLpMw1MlAUotygWCRlrpGBohDlBsUiKXONDBSFKDcoFkmZ\na2SgKES5QbFIylwjA0Uhyg2KRVLmGhkoClFuUCySMtfIQFGIcoNikZS5RgaKQpQbFIukzDUy\nUBQgNxDMUO1opFjPIyzScGSgKEBuIBlmqHY0UqznERZpODJQFCA3kAwzVDsaKdbzCIs0HBko\nCpAbSIYZqh2NFOt5hEUajgwUBcgNJMMM1Y5GivU8wiINRwaKAuQGkmGGakcjxXoeYZGGIwNF\nAXIDyTBDtaORYj2PsEjDkYGiALmBZJih2tFIsZ5HWKThyEBRgNxAMsxQ7WikWM8jLNJwZKAo\nQG4gGWaodjRSrOcRFmk4MlAUIDeQDDNUOxop1vMIizQcGSgKkBtIhhmqHY0U63mERRqODBQF\nyA0kwwzVjkaK9TzCIg1HBooC5AaSYYZqRyPFeh5hkYYjA0UBcgPJMEO1o5FiPY+wSMORgaIA\nuYFkmKHa0UixnkdYpOHIQFGA3EAyzFDtaKRYzyMs0nBkoChAbiAZZqh2NFKs5xEWaTgyUBQg\nN5AMM1Q7GinW8wiLNBwZKAqQG0iGGaodjRTreYRFGo4MFAXIDSTDDNWORor1PMIiDUcGigLk\nBpJhhmpHI8V6HmGRhiMDRQFyA8kwQ7WjkWI9j7BIw5GBogC5gWSYodrRSLGeR1ik4chAUYDc\nQDLMUO1opFjPIyzScGSgKEBuIBlmqHY0UqznERZpODJQFCA3kAwzVDsaKdbzCIs0HBkoCpAb\nSIYZqh2NFOt5hEUajgwUBcgNJMMM1Y5GivU8wiINRwaKAuQGkmGGakcjxXoeYZGGIwNFAXID\nyTBDtaORYj2PsEjDkYGiALmBZJih2tFIsZ5HWKThyEBRgNxAMsxQ7WikWM8jLNJwZKAoQG4g\nGWaodjRSrOcRFmk4MlAUIDeQDDNUOxop1vMIizQcGSgKkBtIhhmqHY0U63mERRqODBQFyA0k\nwwzVjkaK9TzCIg1HBooC5AaSYYZqRyPFeh5hkYYjA0UBcgPJMEO1o5FiPY+wSMORgaIAuYFk\nmKHa0UixnkdYpOHIQFGA3EAyzFDtaKRYzyMs0nBkoChAbiAZZqh2NFKs5xEWaTgyUBQgN5AM\nM1Q7GinW8wiLNBwZKAqQG0iGGaodjRTreYRFGo4MFAXIDSTDDNWORor1PMIiDUcGigLkBpJh\nhmpHI8V6HmGRhiMDRQFyA8kwQ7WjkWI9j7BIw5GBogC5gWSYodrRSLGeR1ik4chAUYDcQDLM\nUO1opFjPIyzScGSgKEBuIBlmqHY0UqznERZpODJQFCA3kAwzVDsaKdbzCIs0HBkoCpAbSIYZ\nqh2NFOt5hEUajgwUBcgNJMMM1Y5GivU8wiINRwaKAuQGkmGGakcjxXoeYZGGIwNFAXIDyTBD\ntaORYj2PsEjDkYGiALmBZJih2tFIsZ5HWKThyEBRgNxAMsxQ7WikWM8jLNJwZKAoQG4gGWao\ndvwwI8AAAAZLSURBVDRSrOcRFmk4MlAUIDeQDDNUOxop1vMIizQcGSgKkBtIhhmqHe3mS369\nT5s109nXjUdYpOHI3ByTe1zjWUO1o914ucW46ZjcdoRFGo7MjTG5zzWeNVQ72o2XmzWjz+/N\nWz/zUTO76QiLNByZG2Nyn2s8a6h2tBsvN2q+929/N6ObjrBIw5G5MSb3ucazhmpHu/VyzaV3\n2o8ccFlD5IncmP3zYb7xcj1ukUTq84f7SPOfzVv/eR9JpD4337xNDm4ixwvSksjw+MPPkWab\nnyONpu//8XMkkfo84JkNIvWxSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTA\nIokAWCQRAIskAmCRRAAskgiARRIBsEgiAM8s0pN+CZMMnieG9hJPLdITzz5Dlh3dXCbLzRaL\ntCfLjm4uk+Vmi0Xak2VHN5fJcrPFIu3JsqOby2S52WKR9mTZ0c1lstxssUh7suzo5jJZbrZY\npD1ZdnRzmSw3WyzSniw7urlMlpstFmlPlh3dXCbLzRaLtCfLjm4uk+Vmi0Xak2VHN5fJcrPF\nIu3JsqOby2S52ZLoSWRwWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFE\nACySCIBFEgGwSCIAFkkE4KlF+hg3o9nimQ5O+Aj5d2U2ylpMyl7WxIVmyzM3NNv8ZYFRzlK+\nQ/7OwWSzmPGzbexJ2cuauNC0PHFD383bYv2P3dvzLBzzPcoIzFcz+l6b+Xq2kZaUvayJC82O\nJ25ouj075kr6aCYZXmbNfPXys3l/tpEtMXtZkxaaPc83FLOTZhbiZdr8LNf/9k6fbWRLzF4O\niDP0/CItmsmzLbR8p1w/Tda/ujF76cgJzZ6nb+hj831MCBmBCSvSMsrKhqjQbHn2hn5GId/A\nbMgIjEX6D7JCs+XJG1qMom6jMwJjkf5NWGi2PGFDh3+YevL8H5Yc2skIzMgi/ZOA0PzmqUX6\nGU9+Hn/+CXlF2j5q95PyqN0yZS9bIkLzm2duaB732EtGYN43d6XnzezZRvZk7GVDXmi2PHFD\nP3kryQhM2jMbUvayJjA0W564obemOfy2KoEQL+PNWoISE7KXZWRotjzRUJO3kxAvi82zv5/t\n4oCQvSwjQ7MlzpDIELFIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBY\nJBEAiyQCYJFEACySCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTA\nIokAWCQRAIskAmCRRAAskgiARRIBsEjhHPxtuot/pu4//n7d/PovlVtxr+H8vUjjw89apDvh\nXofDjUVqLNIDcK/DwSIF414DaJrFuJmu3vgYN6OPzYfmk6aZzJe75M9GzWzz5vbdzcv5tGn/\n9vn23d0lWnbv7/8G+F5D7oB7DaBpVp1YVWK6Cf1k9ZGPzVvNR1uSyfqd6XGR3rdfMmvf7S6x\nZf/+rkidhtwB9xrAqjyL1av5+tVi0qxuRkbN93L52Yy3JflsRt/L79FxkZrmc/2Zpn23u8SW\nY4UjDbkD7jWApvlav5o26zot1t/kNc18/7n1J9afn598a7f/dNur+YnmocKRhtwB9xpAm+5m\nx+r+zOq7sO/v3ef2nz8u0s/8fbIvUneJLccKRxpyB9xrAL+LtHxffRPWjH7+VaTJ/ouPL9Fy\npGCR7o17DeBsyOez8e4eztkivTXjj/lPV6T9Jc4pWKR7414DaNM9Pbmbs2/O9hNf+yJ1bx0V\nafm7J133DjTkDrjXANp0bx5YW36sH2wYbx+Sa29P5t0jbuPmY/3AXrN9hOK7u4/UXWLLocL6\n+725j9rdFfcawC7d23s96/s1n9v7P1/t5zY/YHrb/7ho2j66cPg13SW2F+neH68kjzTkDrjX\nAPbp/liF/m3zeMHmeQn7VqwfOZjtH1N42771tv6K+fbB8sNLtBfZv/813hTpQEN43KsIgEUS\nAbBIIgAWSQTAIokAWCQRAIskAmCRRAAskgiARRIBsEgiABZJBMAiiQBYJBEAiyQCYJFEACyS\nCIBFEgGwSCIAFkkEwCKJAFgkEQCLJAJgkUQALJIIgEUSAbBIIgAWSQTAIokAWCQRAIskAvB/\nKwTeyulMOTIAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"Histogram of residuals.std\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Obtain the residuals from the ANOVA fit\n",
    "residuals=anova.fit$residuals;\n",
    "hist(residuals)\n",
    "# Semistudentized residuals\n",
    "residuals.semistd=anova.fit$residuals/sqrt(mse);\n",
    "hist(residuals.semistd)\n",
    "# Studentized residuals \n",
    "weights=1-1/ns[as.numeric(Spock$Judge)];\n",
    "residuals.std=anova.fit$residuals/sqrt(mse)/sqrt(weights);\n",
    "hist(residuals.std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With residuals, we usually begin the diagnostics with visualizations to explore how the residuals are related with other vairbales.\n",
    "\n",
    "1. $e^*$ or $r$ v.s. $\\hat{Y}_{ij}$\n",
    "2. $e^*$ or $r$ v.s. indices or other structures \n",
    "3. Quantile-Quantile plot or stem-leaf plot \n",
    "4. $e^*$ or $r$ v.s. missing variables.\n",
    "5. ...\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  The decimal point is at the |\n",
      "\n",
      "  -10 | 3\n",
      "   -8 | 21\n",
      "   -6 | 61\n",
      "   -4 | 977\n",
      "   -2 | 3\n",
      "   -0 | 66630941\n",
      "    0 | 4649\n",
      "    2 | 78014\n",
      "    4 | 0771\n",
      "    6 | \n",
      "    8 | 54\n",
      "   10 | \n",
      "   12 | 0\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAWiklEQVR4nO3di1biyAKG0eIiIiJ5/7cdE2y1LwN050+oInuvs+YwM+2krPhJ\nUglQOmC0cu8BwCMQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\nICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJ\nAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAh\nQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkC\nhAQBQoIAIUGAkCBASBAgJAiYIaQCjfmHn/J8OHfYBCQJCQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQICS40aWXHQkJbnL5FXxCgltceS2skOAG115VLiS4gZBgvKtv\ndCIkuE5IECAkCBASJFhsgAAhQYILspDgFiGIcNMqTExIECAkCBASBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQMGtIr8/b4VM4t7vXqTYBdzFjSKf1t09Y30yyCbiTGUPaldXL\ncXj0dliV3RSbgDuZMaRVOX4+PpbVFJuAO5kxpJ8+Wv3/P2d9zCbgTjwjQcC850iHt+GRcyQe\nzZzL35tvq3br0ySbgPuY9zrSbriOtNo+u47EY3FnAwTUE1L5bppNwFTqCWnmTUCSkCBASBAw\n650NN58GCYnGzBjSXkg8rDkP7Y6ryy+eCGwC7mPWc6Tj5RuDEpuAu5h3sWH/7b7ViTbBXSz+\n4p9VO8ZzHV1IjOeOFCExnnu7OiExnpA6ITHa7ZcHH5mQGElIPSExkpB6QmIkIfWExFg66oTE\neELqhESAjoREgo6ERMTCMxISRAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkCmgpp8fehUK2G\nQnJnJPVqJyT36lOxZkLy6jFqJiQIaCUk77BB1YQEAUKCACFBQCshWWygakKCgGZCckGWmrUT\nkluEqFhDIblplXo1FRLUSkgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFC\nggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBA\nSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIE\nCAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKC\nACFBgJAgQEgklHf3HsNdCYnxyod7j+OOhMRopShJSIxVipKExGhC6oTEaKUo6R4h7delbA+T\nboI5Cak3Y0jnSd6c53s3ySa4AyH15g5pV3anrnvblf0Um+AOhNSbO6RVOfWPT2U9xSa4Bx11\n84f0Y6p/n3K/11oV33Et/gTMHdLTj5BWU2yCu8h21OYv01lD2j7vD+Xl/eFpd3m1obFJXLwp\nOmqspFlD+pyfUlanKTbBvUQP61osac7rSMfjfr/dDksOu4sdCWm5hJTU1BQS1OxaupCoiZCi\nWppBkoQU1dIMkiSkqJZmkKhGOxISdRFSUlNTSFSbHQmJ2jTZkZCoT3sZCQkihAQBQoIAIUGA\nkCCgqZBaXM1hGRoKqc3rCyxDOyE1esWbZWgmpFbvwWIZhAQBrYTU7OtUWAYhQYCQIEBIENBK\nSBYbqJqQIKCZkFyQpWbthOQWISrWUEhuWqVeTYUEtRIS1WnxyENIVKbNc2EhUZdGV2eFRFVa\nvV4oJKoipKSmppCgZu+pFBI1EVJUSzNIkpCi3NmwVEKKcq/dYjXaUUMhtTm/SxHbL0JK8nqk\ntiR3TKO7WUiMlt0zbe7lVkJq9iR0AeJ7psVdLCTGsmc6ITFafte0uIuFxEjpXdPmTm4lJMcP\n1QqH1OhuFhIjZUNqdT83E1Krv6kWYKInpLZ2dDshNXrsvABTPSE1tasbCqnN1ZwFEFLXWEhU\nSUidkBgt+rMvpKiWZnDxsj/7jXYkJMYSUk9IjBQ+GmuzIyExWvhJpMmOhMRo8aOx9jISEgGN\nHo1FjQ1pv+66t3VZv6YG9PsmqJ6OxoZ06Cdv1c9htKQl75E2LTyj0SFtykt3LOvupWxiQ+qE\nRHNGhtT/GjqW3flBjpBoTCCkbTkIiYUbfWh3PJRV59COhRu/2FDKc/+EdIgNqRMSzRm9/L3q\nz5C69UtoPH/YBNTPBVkIEBLc6NLFshEh/XzXr1U7HtvlH3QhwS2u/KQ7tIMbXHvOEBLcYK6Q\nXrd//x/6y03A3Vw9jRkb0s45EgswdUhfHbmzgQc2dUir8tJtytvbxuuReGRTh9T/N5/fn42O\nblrloU282ND/Rw9l72UUPLiJQ9q+H9q9lXX3KiQe27QXZIf3bNj0//Wnfxve9U2wOOlF4IzL\ny9Njl7+f+797KsNrKXIqnEbmMsX1lIxLo3JnA3W5cghVKyFRlWsn9bUSElVZaEgTfddNTSFB\nVy981kpI1GSpIX143URv/hbSYi08pO7kOhIJSw/JLUJkNNpRKqR9/3arOU1NIUkLDenrm36O\nDakT0pK12VEqpPU+NqJfN8HCNNmRC7LUp72MhAQR3iASAoQEAWMP7bar/u2DXlfR67FCojWj\n347rOPz/MfvKPiHRmMCbn/z8IEJINGb0+9r9eEZyZwNLNvrQbtW/M+RhNcudDS1eX2AZxi42\nbD7W7GZ4GUWbV7xZhtEXZF+2fUbRd/7+86gavQeLZWjmzoZW7wpmGYQEAaPubJjxPRumu4sC\nAoQEAa0c2gmJqgkJAsaGtF933du6rKMf2GexgdaMDGn4WJdV/8M98UdfekaiaiND2pSX7ljW\n3cvUH30pJKoWuPt7eAmFVTsWLRDSthyExMKNPrQ7HvpXUEx+aGexgaqNX2wY3huylOhtq0Ki\nMaOXv1fDi8zXL6Hx/GETP/6ZjqhXKxdkO69HomazhvT6vB1S2O6uXHXyCtnWLH7XjA7psB1W\n7t6uf91p/e005/LaxLL3SXscLGReav7+z1bXS9qV1cv5rVLeDqvLb9+15D3SIKevo0Pal82p\nn779DZ/Y9+Mdh3pX3nVowTukQRZUu8DbcZ3O12JvmMKf/sjvf94V11bZcV3kzoZbQ/KM9Jjc\ndNIbGdL64xmpv3H1mvdzpMP5TMo50iMRUi9zjvQexg0f2bf5Nt/rU3hU3IuQemNX7bY3LWd/\neN0Nf3y1ff6360jUSEi9yHWkss3eISSkpuioa+oWISolpE5IjCekLhfSMfou+svdHw1yjtQb\nE9LrppTNcG3ouPVBY4slpN6IkF7PM3fs3vr1Bh99uVRC6o0IadPHsyub/lWy24uXheYYFfci\npN6IkM6TVsqqbI8X/vi/WOzuaJKOukhI4XdZ/b4JWiCkLhJScDS/boIm6EhIJOhISEQsPKOR\nIU22XLPsfUKlLv2YCwlucvkHval77RZ//MD9XHnKaCgkZ7Tcz7WDr3ZCssbKHT1MSK76cUdX\n1wOEBNc9TEjTrRDCdUKCACFBwqMsNgiJu3qUkCw2cF9XfvyEBDe5/NMnJLjRpZ89IUFAKyFZ\nbKBqQoIAIVGdFnexkKhMvTvZYgPtqHY3W/6mIdXu5yvDaiaken9TkVRrSNfG1U5IFR87E1Pt\nufADhdTmag5/pdaQro6rqZB4eEKKqmgGmZWQoiqaQeZVZ0dCojGVhvRYiw0sQJ0dCYnWVNnR\nA12QZTEqzKh7mFuE4N4uBS4kCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIA\nIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACGRUOcnscxI\nSIxX6WeDzUlIjFbpp1XOSkiMVevnJ89KSIwlpE5IjFaKkoTEaELqCYmRhNQTEiMJqSckxtJR\nJyTGE1InJAJ0JCQSdCQkIhaekZAgQkgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACHBjS69eZ+Q4CaX305WSHCLK29wLiS4\nwbWP3BAS3EBIMN7VjyWcP6T9upTtYdJNQFhNIZ23vzkPZTfJJmAa1YW0K7tT173tyn6KTcA0\nqgtpVU7941NZT7EJmEhFiw3DAH6M4vLnuwmJytQW0tOPkFZTbAKmUs8F2VK2z/tDeXl/eNpd\nXm0QEtWp5hahbwMpZXWaYhMwoVpuWj0e9/vtdlhy2P3e0bVlEaiYOxsgQEgQcI+Qrh+5CYnG\nCAkChAQBQoIAIUGAkCDA8jcECAkChER1WrxJTEhUps37LYVEXRq9c1lIVKXV1wAIiaoIKamp\nKSTo6rv11KqpkBqbW/6ekKL+/D4t7c0uf0tIUX98B74Wp5e/JKSoP2yi0fnlL7W6n4VEQmy3\ntLqfWwmp2af8RUjul0b3spAYLbtj2tzJQmKs+J5pcRcLibHsma6dkOytavkd1xMSIwmp10xI\nra7mPD4h9doJqdHVnMcnpF5DIbW5mrMAOuoaC4kqCakTEgE6EhIJOhISEYvI6NI3KSS4yeWn\nXSHBLa6cCAoJbnBtaVJIcAMhwXhXb98QElwnJAgQEgQICRIsNkCAkCDBBVlIcIsQRLhpFSYm\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkC\nhAQBQoIAIUGAkCBASBAgJKpz6R1NayUkKnP5PbZrJSTqcuVTH2rVVEjNzS5/7drnENWqoZBa\nnF7+lpCS/vghnU3OL3/n6me11qqZkBqdX/6OkKKEtFRCivp9E81OMH+l2f0sJKrS6m4WElVp\ndTcLibo0updbCanZ31T8rTb3spCoTov7uJmQWn3KZxnaCanRp3yWoaGQ2nzKZxmaCglqJSQI\nEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQE\nAUKCACFBgJAgQEgQMGtIr8/b4U2Ht7vXqTYBdzFjSKf1tw+U2EyyCbiTGUPaldXLcXj0dliV\n3RSbgDuZMaRVOX4+PpbVFJuAO5kxpJ8+SeLyx0oIicZ4RoKAec+RDm/DI+dIPJo5l78331bt\n1qdJNgH3Me91pN1wHWm1fXYdicfizgYIEBIECAkChAQB9YRUvptmEzCVWe9suLkVIbVm8b/9\nZgxpL6RH5UBi1kO74+ryiycCm+AeHJLPfI50vHxjUGITzM/JbTf3YsP+232rE22C2Qmpq2nV\nbuZNkHL7qe8jExIjCaknJEYSUu8eIV2f7cXujhYJqSckxtJRJyTGE1InJAJ0JCQSdCQkIhae\nkeVviBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBBQaUjQmH/4\nKc+H82+qGcg1rQy0lXG2M9CLqvkuqhnINa0MtJVxtjPQi6r5LqoZyDWtDLSVcbYz0Iuq+S6q\nGcg1rQy0lXG2M9CLqvkuqhnINa0MtJVxtjPQi6r5LqoZyDWtDLSVcbYz0Iuq+S6qGcg1rQy0\nlXG2M9CLqvkuqhnINa0MtJVxtjPQi6r5LqoZyDWtDLSVcbYz0Iuq+S6qGcg1rQy0lXG2M9CL\nqvkuqhnINa0MtJVxtjPQi6r5LqoZyDWtDLSVcbYz0Ise47uAOxMSBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBBQQUj7jzH88/uXz2O/LqvdaXi4W30+rM/XOOue\n0NNTKU/H8+OqJ/Q295/m48euPta933fD4Fb93t4MD9f3HtGffY2z8gldDYMbSqp6Qm9092k+\nrj5D2t55KJccy9Opf/J86rrXsjr2w36995j+5Ns4657QXT/E3TDEqif0VvcOaV82HyHty/Od\nx3LJ9jzIfqy7cnh/9FLncL+Ns+4JXZX+yX3Y91VP6K3uHVLZdZ8h7e88lhv0Y92Wt6723/fn\nkOqf0LLq2pjQq+4d0rH7EdK2HJ7eTznvPJ7LTmXzOd5qTz66j3E2MKG7ofUGJvS6Cgb/GdJg\nc+fRXLTvD0Ia2O/DOKuf0JdShswbmNDrKhj85zy+vP8q3dV8PPK26o8+6t/vP8ZZ+YTut6vh\nvKj+Cb1BBYP/aQJPFa+CnlbDb/fq9/vHOD/+puIJ7bqnPvPqJ/QWFQz+5wmseDo35x/JVe37\nffNTOvWOs+szXzUwobeoYPCNhPS23rwND86LTG+1LjJ9jvNDtRM6+FoGrXZCb1PBLH/s6fOF\nhWqn8/B51v48XPY4lDrXw77GWfeE/hjduvYJvVE9Ie36iTydL87V5+1r9avqC/Hfxln3hA53\nNpy2/TlS1RN6q3pCOp1vvqr019JT+bpzbV3xsvK3cdY9oR/32g2zWPOE3qqekN5/ea7Kuta1\n2vItpNNws/K9R/Rnv46z2gkdbvn+GF3NE3qrCkKC9gkJAoQEAUKCACFBgJAgQEgQICQIEBIE\nCAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKC\nACFBgJAgQEgQICQIEBIECKkG3z5nb/isvfPHvv7Ph79e/pjyuj/E/HGZ9hr8EtJ62Cvr/9k3\nQqqRaa/BLz/957/9vySEVCPTXgMhNc+01+DbT//7w/NB3o/PJu/267L6+Gzy3arsvj4FvqyH\n/1+XU3fYlo/PBf88z+p+/erDppTN/5x3MZaQanAppO3wYNP/u03/aPv5hzfl7f2vb+//7vl8\nhrXrfgvp66v35z+zn/UbWw4h1eBrreG3Z5RD2Zy606a8P5W8lNWxO64+Q3opz+9/fX7/V6W8\n9H/7+9d/++pVOfZ/Zn2P728BhFSDCyFt3w/c+sO4bf/wtevb+NxnQxVfi3t/+PpvX12Kw7oJ\nCakGPx/a/RTSr5eYvv/hp/dju7fhgK57Ozxv/hDSt6/evR8VHo+zfUtLI6Qa/GNIr+/Hdrvh\naWrzP89o3766e34/KCyrt9m+qWURUg0uhvTrn/r2T1br/n/9U9N6f3j7Y0jfN3PYrZ0jTURI\nNbgQ0vbr1Ob88PXbH96V/bDgMPyjX0J6PZ8j/XJi5DLTRMxrDX4L6a378ddhqa7b98sFh19W\n7YZ2htWE8n58d/w6R1qXfb9UV3766vV5Zc8z0jSEVINfQlq/n8v8+OvH6c9wbjNcFHr6/qSy\nPl9g2n2cCb2ev37/db3p66tfPv8IExBSDX4J6XXdJ3T+63BvQnk6rxE8f7+zoffycej2Xtfm\n9XBe5T7/uaevOxs+vnq4s0FHExESBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIA\nIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\n8B/QC3rrx1hF3AAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO3di5qqIBSGYSzHznn/dzupHdSsVBawgO99nr23e8xA5R8FtUwN\nwJoJXQEgBQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQ\nAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQ\nQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBMkxY8zhOWX9VsMfHP8K\nY4q/4/tLp+eYu6I8fq3Rn109s0SQHGva7XPK+q36/71sH8HYXIYv/DTHvFSfa3QqaBTLsc0c\ne7Za6SBdil4uBnn5OKf3Y3P6WCPremaJbebYqzULB6k56uxub3zZNUeees6cxxscb0n7+1gj\ngrQG28yxNkjlfar70fHvNvnovzQx2zTHrGb2fmM2t0PFvjDbUzf7UDZhqN6TeLofVe6Txzlz\nnm9w7qamanQ/YAmtfT7YYo41ORieSD06MGVv/rad6OZcque51/O172di1eOMsZsu58x5vcEw\nSP0aEaR12GKO3RrluTvBujfP8tlNeRynTDuw9/xx8Zq7v4Xr+ozDoH3fGv/5MX1+Dmh8n9M/\nIhUfakSQ1mGLOdY0yltL3T+a7fH2z/5aX3f3sy7TZeUxsW+OUOf2n9vPNoPu1aB92/xn2Eea\nqJHkBsgE28yxplle2wNA10D/2lDV7WHm3pifF3VOg39GbyITpKdz/aVGtiudIbaZY22z3DXd\nlmca2gNQfRn9/95+h//cXnaotuZ3kIrXSdn7nN4yD8fXe0zUSG71s8E2c6xrlrdzqbdmOvr/\nZJAOm0fLr0ctvN8TOg16N+9zenVpbat+dj/WCAuwzRzrmuXh1qDff//3OvzTQTo0Paa//Xi0\nunEfmztvj+1owWu0YmLOsC7j/3+sERZgmzl2b5bbx2Gl/NQjmQrSZnwO9nrfU9fPKZukDc7f\nvsyZDtLHGmEBtplj92Z5fgTp4xjZVJB6C9fjFl62b9S+z+Cw82XOdJCmanQVW/9cECTHHo33\nzwyOTY2//vzJIG3bs7RmtPo9SIM76kx1nTFnOkjjGjVLVzWWIUiOPRrv9RGkZ7v9G8yfDNLp\n/tKiHRMf5eB1j3ejnDHnQ5BGNWoivxVZ95wQJMeejXf/nGofFurdazeYGP5zbno5f+fLY1hu\n+N7Hv0371FFlBrfUfZzzKUjDGjWnhjyRtBRBSsH579Mh5PMciCJIgACCBAggSIAAggQIIEiA\nAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAgg\nSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIE\nCCBIgACCBAjwECQDRGZFK5cPToAiAEkECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQ\nQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQgLm+PFFOkIB52hR9ihJBAuYxvb8/zFzx\nfk4RJKhjRv9Oz13xhg4RJKhDkAABBAmQQB8JEMCoHSCC60iAWwQJEECQAAEECRBAkAABBAkQ\nQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEE\nCRBAkAB7nz+n68sy8tUIUAQg5hYjggRYMjWndoAt8/xrxXJuESRE4+tn6/9e0CmChFiYwT8r\nlnSJICEOz9E6ggSsZiamVizsDEFCDMzk5IqlXSFIiID5ML1icUcIEvQzH/+zYnk3CBK0G90U\nRJCAFcZNlCABy721UK9BOu1K0yirk6siAA/eG6jHIF035mXrpAjAg6lnJjwGqTLF4dxOXY6F\nqVwUAbg32To9Bqkw5+f02RQuigCcm26cHoM0OCB+f6KQIEGrD22TIxKwwKem6bePdLy0U/SR\nEKfPJ1I+h7+3vVG7zdVJEYBDX9ql3+tIVXsdqSh3XEdCfL41S+5sAOaRHh9z1MpNn5sigPW+\nN0o9QfJcBLDEr1/uBAn47WeLJEjAT78bpNc7G2Z3gwgSNJnRHj0GaU+QEKU5zdHnqd25+P7w\nhEARgLR5Y8he+0jn7zcGSRQByJrZFv0ONux79606KgKQNLcpMmoHfDa7JRIk4JMFt9gQJORg\n1W1nSxYhSEhfm6LFUVr0eoKE9Jne3wsXcvPqtYsoLAL5MKN/Fy3k5uUrF1FYBPKxIkhuzwNX\nL6KwCORjeZCWN0CChPQt7SP5SQVBQmQWjtqtaX4ECTlYch1pVesjSEDfyg8MIUhAz+rHVv0V\npasIYMLqlkeQgKf1DY8gAXc2n6dIkICOVasjSEDLrtERJKBh2eYIElDbNzmCBFgNM9zfwcsi\nCosAngTaG0FC9iSaG0FC7kRaG0FC5mQaG0FC1qS+HZIgIW52X5Qq1tIIEmK27hPrXovLVSRk\n4SGLQBJWfWLdaGnBioQqPWARSMGqT6x7W1oCQULEbIIkNcywugYECVpYBEm4jYWvQaAikITV\nfSTpJqagCmGKQBLWjtqJtzCChLituY4k2z3q3tLLIgqLQL5cNC+ChNw4aV0ECZlx07gIEvLi\nqG0RJOTEwTDD/Y29LKKwCOTIXcMiSMiHw3ZFkJANl82KICEXTlsVQUIenA0z3N/eyyIKi0Be\nXDcpgoQcOG9RBAkZcN+gCBKS57h71JXhZRGFRSAbXloTQULi/DQmgoS0eWpLBAlJ89WUCBIS\n5mOY4V6Sl0UUFoEMeGxHBAnJ8tmMCBJS5bUVESQkym8jIkhIkr9hhnt5XhZRWASS5r0FESQk\nyH8DIkhIT4D2Q5CQnBDNhyAhMb6HGe6lellEYRFIVKC2Q5CQlFBNhyAhJcFaDkFCOsJ0j7qi\nvSyisAikJ2SzIUhIRdBWQ5CQiLCNhiAhDYHbDEFCCgIOM9wr4GURhUUgJeEbDEFC/BS0F4KE\n6GloLgQJsVPRWggS4hZ8mKFDkBA1LU2FICFmaloKQULE9DQUgoRoKeketQgSYqWqlRAkREpX\nIyFIiJOyNkKQECVtTYQgIUKahhk6BAnxUdg+CBKio7F5ECTERmXrIEiIjM7GQZAQFX3DDB2C\nhJiobRkECRHR2zAIEuKhuF0QJERDc7MgSIiE1mGGDkFCHJS3CYKEKGhvEgQJMVDfIggS9NPd\nPWoRJKgXQ3MgSNAuitZAkKBcHI2BIEG3SNoCQYJmEQwzdAgSFIunIRAk6BVROyBIUCumZkCQ\noFVUrYAgQadohhk6BAkqxdYECBI0iq4FECQoNG4ARv2JHkGCOuPUtP9XHiWCBG3e9r6Z/rEq\n/oO03xhTHp0WgZh9yJHyVuExSN2xeWtalZMiEL/3fU+QRss1C1amutb1pTJ7F0UgehO7niCN\nlmsWLMy1mb6ajYsiELnpEQX6SMPl+mMv71vM9K0sAnH7sN8ZtRsu1yz49whS4aIIRO3zbtf/\nu9VrkMrd/mgOt8lr9X20QftWgwtR73WvQXqethlTXF0UgagMjzNm4mfx8Hkd6Xze78uyHXKo\nvuaIIOVg2PPpJmLoDU3jzgYEMhiLMxM/iwpBQhiDq0Nm4mdxIUgIox+acYAi3P8ECWH0QmMm\nfhYbgoRAHv2h/tgCfSRZMW5JLHQfoTMTPwtTISsECcE0KYrvWdhpBAkhJbOnCRLCifToM4Ug\nIZiUdjNBQihJ7WWChEDS2skECWEkto8JEkJIaJihQ5AQQHo7mCDBvwT3L0GCdynuXoIEXx53\n/yS5dwkS/Hjcj5rcMEOHIMGPx2ebBK6GKwQJXjweJk911xIkePF8jC9wPVwhSPDi8cReqruW\nIMEPM/iYk+TYBmm/qevLxmxOUhV6LwJJiPcp8lksg3RstkzRfBCxaJKS3dwZM9E+RT6LZZC2\n5lCfzaY+mK1YlWqClKDUd6llkJrfMefmmyVkf9mkvtWzk/KxqCMQpNIcCRK+yWB/Wp/anY/N\nd4ZxaofPctid9oMNxuyaA9JRrEp1Hls+H1nsTevh76L97r3NQag+E0UgbnnsTC7IwhmT7s3e\nbwgSHEn8CuyIdZCOZTtydxGqz1QRiFJ3R1Auu9I2SNvu65VNIZqkXLZ+ykzvTwYsg7Q322sT\npL35E6tSnc3GT9qje5TJvrQMUmGuLs6FM9n4SXu0iEz2pcCdDQQJE1L+pJMJlkHa3I9IzY2r\ngnLZ+gkzjNotWeTeRzoWZi9WpZogRcsMTuiSfm5ixHbUrjQd0VvtCFKcnsegjAL0IHIdyZSy\ndwgRpHj0jzqpf8DJN9zZAAuDftBjr2W59wgSLAxG5h7/ZLn3LIJkhgLXCgEMj0E5n9kRJFgY\nncx1wwx57jxO7bDU67fmOEiZjtg1CBKWmRhfMK95mcbI8tRucHoXuFbwZDi+MLh7IefdRpAw\nh/k0xt3b81nvNU7t8FvvwPPx6JPtSV2HIOG33uncpyDlvssEHqNoFYVEbaaKQHCD8Ew/HZH9\nHhMK0oU+UsKGQZp6OoIdZhGk4+B6LM8jpWvqctFgPvvL6oi06eeIr3VJ2PTp3HBu5qT6SLLY\nNbp8fdiVndVg1A5zfL5QyL5qESRYYVd1bIO023BnQ8YYZniwDNKOW4Ryxn56sgyS8KcHTRUB\ntdhNL4zaYS32Uo9lkEpzFavKhyKgFDupzzJIl2IreiV2ogioxDDDkPWpHYMNWWIPjRAkrMAO\nGuOCLJZj/7whSFiK7tEEqSCdStua/CwCOrBzptgGqaKPlKzpncq+mWQZpFeOjmJVqtlZGnx4\ndIJdM836FqFDvTWXy9bwYF9i+GyGRQRuEdrdjkZn2W8aY3cFN/lpQQwzfCQQpGNz4yp9pMRM\nBYnd8pn1vXaH+mI29YkgJWYiSOyVLyyDdGwCtG0GG/7EqlSzyzR46yOxU76xfkK2+d+fMZVQ\nfSaKQBDjUTv2yVfc2YBPBt+zzC75jiBhBnbILwQJv7E/fuIxCvzE7viNIKFvakeyN2aQObU7\nbUVv/mbXBTJ1fx3DDLMI9ZGuXEdKwcT9deyKeaQGGzi1SwB3M6wnFKS94Rv74vceJHbEXGKD\nDTuxKtXsv0DGQaJ7NJ9QkDayn1zMDgxj2EdiLyzABVm8DEbt2AlLECT0va4jsQ8WsQvSddc8\nQlE2n9dQSnaS2InBsQuWsQrS/tlFuuxFP/2EvRgYwwxL2QTplqOq+cyTc2UKPrMhJWz/xSyC\ndH19BteBj+NKCZt/OYsgVa/HYkvhR2TZkyGx9VewCNLGXO5TF7MzG6ka9YtAAGz8NSyC9OqQ\nnkvutUsFwwzriARp/B9r7MxQ2PIriZza1e1n2wlidwbChl9LZrBhMC2A/RkG2301meHvo5H9\ndnN2aAh0jyzYXJA9GFOd6/aCrOxlJIIUAhvdhtUtQofifotQIZsj9qlf7Z2qbHMrljet7svm\nptW96HldzU716sM3imERHqPIXruxCZIlgpS7x/GIbW6FIOXukSG2uRWClKvHs7Bm6tPssBRB\nylNvhIEcSSBIeXqlxzBqJ4EgZel14/ErS7BBkLLERwVJsw7SsWx+nZWXDy9eh/3r2OOj69jQ\nUmyDtO2+GskUokli/7rGVVhhlkHam+212R17vtYlLu1vP4IkxzJIhbm6uFeLHeweMRJlGaTu\nF1tNkKLDJpZlGaTN/Yh0nvWo+WlXtk9dlO3nSsrWCouwhYXJ9JGOhfn9vS7XTe+rm79/MCu7\n2TE2sDTbUbtyVjBalSkO53bqcgve1894YD87Re9Insh1JFMeZixXmPNz+vz9qzLZ0W7cb2QI\nXY0Uebyz4cfn4Jm+lUXgG24HcshjkDgiBcZt3g5ZBMmYZQeRWx/p2N3/QB8phOHjR5DlMUjd\n7USPryb7+nkp7GsHeKbcJa93f5+qdpCvKHdcR/KPZ8pd4jGKfBhy5I51kA7NQeaPD4jU63na\nzaOwDok8RtFcSZKq0HsRsPFKD/d7u2QZpKr7tOJZtwitLAJWnmN1bFOnrB+j6K4NzbtpdVUR\nsMEz5Z4IPEYxnBDBbpfCM+WeWJ/aPY5Iop0kdrsUrsJ6YjvYsGv7SKdixt3fa4uAje7JSzao\na9andk7uNWW/izHPh5jhEkFKHTfTe8GdDYljU/pBkNLGlvSEIKWMkzpvbINUFdL9o7cisBab\n0R/r60hOng6nBUhgK3pkPWoneo/dVBFYg7E6z6RuEZJFG7DDJ3t7Z31q9/WR8bVoAqu8njzi\nmXLfrJ9H2sp+M9JEEZin/+QRHxjkm22Qjgw2aPEKz+vABF8sg7Rj1E6L3pNHPITkn/WDfYza\nKfFKD2d2ATBql4r+I5Z8zIl31qd2jNppMTgMMfjtmfWDfdsfn/W4Co1gBQ5DIck9jyRWpZog\nLdDf8hyGwiFIURsehchRODxGEbVhvyhgRbJHkGI2uGDERgtJKkgnPo4rgH6Q2GZBWT/YRx8p\noIm7GRCG2IN9ot9HQauY6dFHYpghNOtbhA711lwuWyN6OYlm8WZ6mJvvV9ZC4Bah3e1odDai\nH7VKw2j0AtIf5h4PefM5qhoIBOnY3LhKH0naIC/mw9/jSYRiGaTydmp3MZv6RJCkTaXGjKaH\nEwjIMkjHJkDtt/b9iVWppm3UHy4RTQSJ/pEO1jetNv/7M6YSqs9EEXmaGSS2lBLc2aDUVF6+\n9pQQFEHSavR4UV1PjNqxndSwCdK1aidPG1MIP3BOA3l7vGjiOhLdI0VsglS0e7L7HCG+sU/e\njwutbCRNLIK0N9vmOfOiONfXrTkErlV22EaqWARpa5rPhjyZXfs3dzb4xSbSxSJI3YlH1d1l\nxwVZv9hCylgHadMbShJDM/mOYQZ1LIK0aU7tLt0tDVdTCFaKIH3H5tHHIkhVM9jw1z2ItOcW\nIX/YOgpZBOlaPMe998acBStFU/mGjaOR1QXZxz12hnvt/GHbqCRyi5AphT9ulcbyCcMMSnGv\nXVTYMFoRpJiwXdQiSBFhs+hFkOLBVlGMIMWCYQbVCFIk2CS6EaQ4sEWUI0hRYINoR5AiQPdI\nP4KkH1sjAgRJPTZGDAiSdmyLKBAk5dgUcSBIqjHMEAuCpBnbIRoESTE2QzwIkl5shYgQJLXY\nCDEhSEoxzBAXgqQTWyAyBEml7DdAdAiSRrmvf4QIkj50jyJEkNTJeuWjRZC0yXndI0aQNDE/\nvu0SahEkPZoQkaRIESQ9TBel0NXAGgRJDXNf7SxXPnoESY3HSV2WKx89gqSGGf2LmBAkLTiz\nixpB0qEbsau5rSFWBEmFR/eIGMWKIGmQ2/omiCApkNnqJokghZfX2iaKIIVGtygJBCmwjFY1\naQQprHzWNHEEKahsVjR5BCkgukfpIEjh5LGWmSBIwWSxktkgSKHksI4ZIUiBZLCKWSFIQTDM\nkBqCFELq65chghRA4quXJYLkX9prlymC5F3SK5ctguQZwwxpIkh+pbtmmSNIXiW7YtkjSD6l\nul4gSB7RPUoYQfImyZXCHUHyJcV1whNB8iTBVUIPQfIjvTXCAEHygWGG5BEkDxJbHUwgSO59\nWxs+Nj8RBMm5LyvDF7kkgyC59vV49PMViIT/IO03xpRHp0Uo8vVwY0b/Il4eg9S1qa1pVU6K\nUOf7ihCkdPgOUmWqa11fKrN3UYQ2P9aDIKXDd5AKc22mr2bjoghlfq4GfaRk+A7So8/wfagq\niaY1YzSOUbtk+A7S3yNIhYsiNJm3DlxHSoTXIJW7/dEcbpPX6vtoQwKNK4FVwAJeg9RpJ4ur\niyL0iH8NsIjP60jn835flu2QQ/WeI9O3tggtol8BLMSdDQ7E/4sASxEkeXHXHquECNLvX9hR\nN8WoK4+VCJK0mOuO1QiSsIirDgsESRTDDLkiSJJirTesESRBkVYbAhj+lhNnrSGCIImJstIQ\nQpCEMMyQN4IkI74aQxRBEhFdhSGMIElYUN8Ebm3HBIJkb0E0eLY8VQTJ2pLK8mknqSJItpbn\nKK71wywEydKiqhKkZBEkO8tqSpCSRZBszLhrcPgS+kipIkgWZsSoHqaNUbtUEaT1ftdy4gDE\ndaQ0EaTV5uYoktWBFYK01nsd3w42BCkfBGmd9xO0ie4PQcoHQVplooJTA3IM0mWDIK3xMUdv\nIwv11MEL6SFIK0xV78NpHIN0mSBIi01Hg/5Q3gjSUp/qRn8oawRpoY9Voz+UNYK0zLea0R/K\nGEFaRG3FEBhBWoAjDj4hSPPprBVUIEizqawUlCBIc/XrxLACRgjSTDydh28I0iyD2HDpFW8I\n0hw8HIEfCNIMU89GqKskgiJIv/HcK34iSD+9Pws7/WPkjCD9MDE6x6gd3hCk76ZrwnUkjBCk\nr9RUBMoRpG+01APqEaTPOH/DbATpIxWVQCQI0ica6oBoECS9VUBECJLWGiAqBGmy/NAVQGwI\nkr7iESGCpK10RIkg6SockSJImspGtAjSqGRyhDUIko6CETmCpKFcRI8ghS8WCSBIr0LJEVYj\nSCHLRDIIUrgikRCCFKpEJIUghSkQiSFINcMMsEeQOBxBAEEiRxBAkMgRBGQfJHIECZkHiWEG\nyMg7SMQIQrIOEjmClJyDRI4gJuMg/SqE727BfNkG6VdK+DYxLJFrkH6WYDxVBGnIMkgzTtrM\n6F/gqwyDdEuR+X1m56MmSEeWQTK/iyBIWCS/ID2ORbOSRI4wT4ZBmlcGo3ZYIrsgmdknbVxH\nwny5Bclw0gYX8grS63yNow1EZRWk51kdMYKwnIJEeuBMRkEiR3AnnyCRIziUS5DoFcGpTIJE\njOBWHkEiR3AsiyCRI7iWQZDoHsG99INEjOBB8kEiR/Ah9SCRI3iReJDIEfxIOkgMM8CXlINE\njOBNwkEiR/An3SCRI3iUbJDIEXxKNEgMM8CvNINEjOBZkkEiR/AtxSCRI3iXYJDIEfxLLkgM\nMyCE1IJEjBBEYkEiRwgjrSCRIwSSUpDoHiEYr0E67UrTKKuTgyKIEcLxGKTrxrxsxYsgRwjI\nY5AqUxzO7dTlWJhKuAhyhJA8NtnCnJ/TZ1PIFkGOEJTHIA3GAr4PDCwtgmEGBJbEEYkYITS/\nfaTjpZ0S7iORIwTns1u/7Y3aba5iRZAjhOf3OlLVXkcqyp3gdSRyBAV0DjTPL4JhBqgQeZCI\nEXSIO0jkCEpEHSRyBC30BMn0zVvATUWA5bze2TA7K3OKIEZQxGOQ9qJBIkfQxOep3bn4/vDE\nkiLIEVTx2kc6f78x6HcRzwMZOYIufgcb9r37VpcX0aao+YthBmijZ9TudxHm/jcxgjoRBcl8\nmwkEFV+QOK+DQiGC9DsKX4LEmR00iihIdTfMQI6gUFRBYsQOWsUUpHruXXiAbzEFiRBBrYiC\nRI6gVzTD35zUQbNYgkSMoFokQSJH0C2OIJEjKBdDkOgeQb0IgkSMoJ/+IJEjREB9kMgRYqA9\nSOQIUdAdJIYZEAnVQSJGiIXmIJEjRENxkMgR4qE3SOQIEdEaJIYZEBWlQSJGiIvWIAFxWdHK\n5YNjT2OlFNZJYZXyrZPGFVdZKYV1UlilfOukccVVVkphnRRWKd86aVxxlZVSWCeFVcq3ThpX\nXGWlFNZJYZXyrZPGFVdZKYV1UlilfOukccVVVkphnRRWKd86aVxxlZVSWCeFVcq3ThpXXGWl\nFNZJYZXyrZPGFVdZKYV1UlilfOukccVVVkphnRRWKd86aVxxlZVSWCeFVcq3ThpXHIgOQQIE\nECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBOgM0n5j\niuoauhZje00bqyrYRHP4akrqVrxRtd8IUChrJmdN30+zbTfRJnQ1RlRtopa3pqRtxRtn83dt\nfrn9ha7IwLlQ1EpOpjg3NTqFrsiAqk3U8teUlK14q7x/fbOquu3NVlGFKnO8/X0wu9AV6dO1\niVr+mpKyFe/TtVNMpalCpbnUze/bMnRF+nRtor6sg3Q129BV6DurSrbReNDWtYl6fDQllSve\n2rdnL5ooaiUqg1QrrFDLR1NSueKNS6HqtKWhqJUQpAW8NCWNK964FqpO7FqKWglBms9PU9K0\n4v0vlN4quUTSr5OiVlIQpNn8NCVNK/5qtJfN9hK6Nh2lQepG7S66Ru1qVZvozldTUrfijaOu\nAbsHRa1k13afj6YKXZERRZuo460paVvxxkVnjjS1Ep13NqjaRC1/TUnZirf+jOmfUamhqUKb\ndgOp+4WjaRM1/DUlZSveMgTpp2t793foWrzRtIka/pqSshUH4kSQAAEECRBAkAABBAkQQJAA\nAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBA\nkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkJwyPWu+zu7YvcnH\nd+/951ptjNlU12Xvq+0r9qLFdnTKMkgb073Jx3d/TR8exewXvS9BEsJ2dO/RWJc32h9L9Gbf\nclRd6vpSzUrSa0GCJITt6J6HIF0Lc+ymjsb8PrsjSOLYju71glSZYtdO7zemuB86bpObfTf/\nujFlb+b9+7i75avCbC/NxLE09y80f6Vgb57fcF6Z3WNO+3f/5ZeyLb/3vt0LX7U5bo3ZHh1t\niJQRJPdeQSofnZh2wmybn26fk+38qjezH6T2ZcXtYLPrukJV3Q9Sac6PyVO35KPgwcuLZnL3\nFqRXbfbz+1kYIEjuvYK0vd6a6qY5/7pNXbfN6djBFOf6XJjDff5g5isPh+aHf00eTPPSQztn\n+gztdaDppnovf5Tfnz0osGgCeWhegmUIknuvIJ3u/yvbfsy1OY8r287NsTuMNPP7M18Nvmzm\nXU0xeM85QRq8/PQ2e1QbYzitW4cguTccbOia72NI/DFzMDmaOR6mvhx32wVBGr38/X17BbBa\nXcQAAAG5SURBVFa308vzucZiBMk94SBtH3On+0jn4aHs7eVfg1Tvmm5UcXG0JRJGkNybCtJ4\n5ihT9ecl/sxmf7yMgnQftTtfmoPKoHP19vKpIPVre6w29JGWI0juvQepfHVFHn2k52GkP/O1\nxPbZR2p/Ng7S/TpSacpupKCdc3qm5HuQynHHyNAqFmOTufcepHao7nYYKUejdu3LejONudyX\n2Dcja1U3aneqz+M+UnMdtrmzYXc/L9uYfTMOZyZe3v3wUr+C1Ctw043xcURajCC59x6ke8el\nbfP960jd614zN+Z5DHpeR6ruHZrT8MhxfPZ0motA7fWg8j58MHx58/frfUe1OTxfjGUIknsT\nQWruJTB/XZ9+XzzvbLgv8Jx52jyD1A6otQv83WJ36p8Ldu53fx+37c0Ru8L8dfPHL2/+fr3v\nuDbtnQ3kaDmClJzjLnQNckSQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAA\nAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRBA\nkAABBAkQQJAAAQQJEPAPi7q+6Qci4BMAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"Normal Q-Q Plot\""
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the residuals (or the other two versions) against fitted values\n",
    "plot(residuals~anova.fit$fitted.values,type='p',pch=16,cex=1.5,xlab=\"Fitted values\",ylab=\"Residuals\")\n",
    "\n",
    "# Plot the residual against certain orders\n",
    "# No clear orders make sense in the Spock trial data\n",
    "\n",
    "# Stem-leaf plot  (or use histogram, or qq-plot )\n",
    "stem(residuals)\n",
    "qqnorm(residuals);qqline(residuals)\n",
    "\n",
    "# Plot residuals against missing variables\n",
    "# Not applicable on Spock trial data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 4.1.10 Testing equal variance \n",
    "\n",
    "Estimate $\\sigma_1^2,\\sigma_2^2,\\ldots, \\sigma_r^2$ separately as, for $i=1,\\ldots, r$, \n",
    "$$\n",
    "s_i^2 = \\sum_{j=1}^{n_i} \\frac{\\big(Y_{ij}-\\bar{Y}_{i\\cdot}\\big)^2}{n_i-1}.\n",
    "$$\n",
    "We want to test the null hypothesis $H_0: \\sigma_1=\\cdots =\\sigma_r$ against the alternative $H_a:$ not all $\\sigma$s are equal. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".dl-inline {width: auto; margin:0; padding: 0}\n",
       ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
       ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
       ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
       "</style><dl class=dl-inline><dt>1</dt><dd>25.3894444444444</dd><dt>2</dt><dd>43.3256666666667</dd><dt>3</dt><dd>21.095</dd><dt>4</dt><dd>35.6275</dd></dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[1] 25.3894444444444\n",
       "\\item[2] 43.3256666666667\n",
       "\\item[3] 21.095\n",
       "\\item[4] 35.6275\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "1\n",
       ":   25.38944444444442\n",
       ":   43.32566666666673\n",
       ":   21.0954\n",
       ":   35.6275\n",
       "\n"
      ],
      "text/plain": [
       "       1        2        3        4 \n",
       "25.38944 43.32567 21.09500 35.62750 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate the variances for each group:\n",
    "(vars = tapply(Spock$perc.women,Spock$Judge,var))\n",
    "alpha=0.05;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hartley test.【balance design】** The test statistic is \n",
    "$$\n",
    "H=\\frac{\\max(s_1^2,\\ldots, s_r^2)}{\\min(s_1^2,\\ldots, s_r^2)}.\n",
    "$$\n",
    "At significance level $\\alpha$, reject $H_0$ if $H>H(1-\\alpha;r,n_i-1)$ when all $n_i$s are the same (balanced design). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "8.43999329265814"
      ],
      "text/latex": [
       "8.43999329265814"
      ],
      "text/markdown": [
       "8.43999329265814"
      ],
      "text/plain": [
       "[1] 8.439993"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "7.18532327249131"
      ],
      "text/latex": [
       "7.18532327249131"
      ],
      "text/markdown": [
       "7.18532327249131"
      ],
      "text/plain": [
       "[1] 7.185323"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Hartley test:\n",
    "H.stat=max(vars)/min(vars);\n",
    "#install.packages('SuppDists')\n",
    "library(SuppDists) # The distribution is in this package\n",
    "# Both df and k only take integers:\n",
    "qmaxFratio(1-alpha,df=floor(sum(ns)/length(ns)-1),k=length(ns))\n",
    "qmaxFratio(1-alpha,df=ceiling(sum(ns)/length(ns)-1),k=length(ns))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bartlett test.【large sample sizes in each cell】** The test statistics is \n",
    "$$\n",
    "K^2=(n_T-r) \\log({\\rm MSE})-\\sum_{i=1}^r (n_i-1) \\log (s_i^2).\n",
    "$$\n",
    "这个很像likelihood ratio test\n",
    "\n",
    "We know that $K^2 \\geq 0$ from Jensen's inequality. Under $H_0$, $K^2$ is approximately $\\chi^2_{r-1}$ assuming that $n_i$ are not small. Reject $H_0$ if $K^2> \\chi^2(1-\\alpha;r-1)$ at significance level $\\alpha$. Related to the likelihood ratio test. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "7.81472790325118"
      ],
      "text/latex": [
       "7.81472790325118"
      ],
      "text/markdown": [
       "7.81472790325118"
      ],
      "text/plain": [
       "[1] 7.814728"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Bartlett test:\n",
    "K.stat= (sum(ns)-length(ns))*log(mse)-sum( (ns-1)*log(vars) );\n",
    "qchisq(1-alpha,df=length(ns)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Levene test.** \n",
    "1. Create new data with $d_{ij}=|Y_{ij}-\\bar{Y}_{i\\cdot}|$. \n",
    "2. Treat $\\{d_{ij}\\}$ as response variables\n",
    "3. Calculate the $F$-statistic for $H_0: \\mathbb{E}[d_{1\\cdot}]=\\mathbb{E}[d_{2\\cdot}] = \\cdots =\\mathbb{E}[d_{r\\cdot}]$\n",
    "Reject $H_0$ if $F^*>F(1-\\alpha; r-1, n_T-r)$ at significance level $\\alpha$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value Pr(>F)\n",
       "Judge        3   5.64    1.88   0.173  0.914\n",
       "Residuals   29 314.70   10.85               "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Levene test:\n",
    "Spock$res.abs=abs(anova.fit$residuals);\n",
    "summary(aov(res.abs~Judge,data=Spock))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1.11 Remedies for departures from model assumptions \n",
    "\n",
    "**Weighted least squares**\n",
    "\n",
    "Idea: $\\sqrt{w_i}\\epsilon_{ij} \\sim N(0,1)$ if $w_i = 1/\\sigma_i^2$. \n",
    "\n",
    "Find weighted least squares estimator by minimizing $\\sum_{i=1}^r \\sum_{j=1}^{n_i} w_i(Y_{ij}-\\mu)^2$ for the common mean \n",
    "$\\tilde{\\mu}=\\sum n_i w_i \\bar{Y}_{i\\cdot}/\\sum n_i w_i$. We will have a new set of ${\\rm SSTR}$ and ${\\rm SSE}$. In the end, we still have $F^* \\sim F(r-1,n_T-r)$. \n",
    "\n",
    "In practice, we plug in $w_i = 1/s_i^2$, and the null distribution remains the same.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Rank test**\n",
    "$$\n",
    "F^*=\\frac{ {\\rm MSTR}(R) }{ {\\rm MSE} (R)} = \\frac{\\sum \\sum (\\bar{R}_{i\\cdot }-\\bar{R})^2 /(r-1)}{\\sum\\sum (R_{ij}-\\bar{R}_{i\\cdot})^2/(n_T-r) } \\sim F(r-1, n_T-r),\n",
    "$$\n",
    "where $R_{ij}$ is the rank of $Y_{ij}$ among all $n_T$ observations. Works when the sample size is large.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "            Df Sum Sq Mean Sq F value   Pr(>F)    \n",
       "Judge        3   1846   615.4    15.6 3.15e-06 ***\n",
       "Residuals   29   1144    39.5                     \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The rank test\n",
    "Spock$rank.perc=rank(Spock$perc.women)\n",
    "summary(aov(rank.perc~Judge,data=Spock))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Kruskal-Wallis test**\n",
    "$$\n",
    "F^*=(n_T-1)\\frac{ \\sum_{i=1}^r (\\bar{R}_{i\\cdot}-\\bar{R}_{\\cdot \\cdot})^2 }{ \\sum \\sum (R_{ij}-\\bar{R}_{\\cdot\\cdot})^2} \\sim \\chi^2_{r-1},\n",
    "$$\n",
    "if $n_T$ is large. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tKruskal-Wallis rank sum test\n",
       "\n",
       "data:  perc.women by Judge\n",
       "Kruskal-Wallis chi-squared = 19.757, df = 3, p-value = 0.0001906\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Krusal-Wallis test:\n",
    "kruskal.test(perc.women~Judge,data=Spock)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Box-Cox transformation**\n",
    "\n",
    "$$\n",
    " Y(\\lambda)=\\frac{Y^{\\lambda}-1}{\\lambda}, \n",
    "$$\n",
    "and $Y(0)\\equiv \\log(Y)$ for $\\lambda=0$. \n",
    "\n",
    "To tune the parameter $\\lambda$, we can calculate the likelihood $L(\\lambda) \\equiv  \\max_{\\mu,\\sigma} L(\\lambda;\\mu,\\sigma)$. Then $\\lambda^* = \\arg \\max L(\\lambda)$. \n",
    "\n",
    "It can be shown that $\\max L(\\lambda)$ is equivalent to $\\min {\\rm SSE}[Y^*(\\lambda)]$, where\n",
    "$$\n",
    "Y^*_{ij}(\\lambda)\\equiv \\begin{cases}\n",
    "\\frac{Y_{ij}^\\lambda-1}{\\lambda \\dot{Y}^\\lambda -1} & \\lambda\\neq 0\\\\\n",
    "\\dot{Y}\\log(Y_{ij}) & \\lambda =0\n",
    "\\end{cases},\n",
    "$$\n",
    "where $\\dot{Y}$ is the geometric mean of $Y$. \n",
    "\n",
    "Box-Cox for equal variance amounts to minimize the Bartlett statistics, or other test statistics. We can use the `boxcox()` function in library `MASS` for Box-Cox transformation in `R`. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1.12 Power calculation\n",
    "\n",
    "\n",
    "Recall that type II error is *accepting the null hypothesis when it is false*. The type II error rate is the probability of type II error given repeated samples from the true population, i.e., ${\\rm Prob}({\\rm Not\\  rejecting} \\ H_0 | H_0 \\ {\\rm false})$. The **power** of a testing procedure is defined as *the probability of rejecting $H_0$ when it is false*. We have the following equality \n",
    "$${\\rm Power} = {\\rm Prob}({\\rm Reject}\\  H_0 \\mid H_0 \\ {\\rm false}) = 1 - {\\rm type \\ II \\ error \\ rate}. $$\n",
    "\n",
    "The power of a test depends on multiple factors including \n",
    "\n",
    "- Deviation of the alternative $H_{a}$ from the null $H_0$, or the strength of signal.\n",
    "- The noise lelve, or the degree of random fluctuation. For instance, a larger ${\\rm var}(\\epsilon)$ in an ANOVA model leads to lower power of the F-test. \n",
    "- Significance level. Recalling that the significance level is the type I error rate. A larger significance level results in larger power, which reflects the trade-off between type I and II error rates.\n",
    "- Sample size. The more samples we have, the more powerful the test is.\n",
    "\n",
    "Among the four factors above, we can only change the significance level during the analysis. However, given the relationship between type I and type II error rates, simply increasing the significance level will decrease the reliability of the test results. As an extreme case, a test at significance level $100\\%$ has a power of $100\\%$, but this test is meaningless. Moreover, changing the significance level during analysis is considered as a form of p-hacking. \n",
    "\n",
    "Prior to the study, we can determine the sample size of the study to achieve a desirable power at a given significance level. This is known as the **power calculation**. To this end, we need to know the other factors that influence the power. Here the significance level is pre-specified and hence known to us prior to the study. Deviation of the alternative and the noise level, however, are population quantities that are unknown, which is the reason why we need to conduct a study at the first place. In practice, we can either conduct a pilot study to collect a few samples to estimate these parameters, or to specific a level of clinical significance. Finally, we will need to use the property of the testing procedure to \"predict\" its power. Briefly, we need to obtain the distribution of the test statistic under the alternative hypothesis, and then to calculate the probability that the test statistic surpasses the critical value. The distribution of F-test in Section 4.1.6 will come in handy if we want to calculate the power before hand. "
   ]
  }
 ],
 "metadata": {
  "author": "\\mu",
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}