{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np # 导入包\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "nsample = 50 #生成50个样本点\n", "x = np.linspace(0, 10, nsample) # 从0-10之间生成50个数\n", "X = np.column_stack((x, x**2)) # 生成第一列是1，第二列是x，第三列是x的次方\n", "X = sm.add_constant(X)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4.72063062, 2.098422 , 2.98932883])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta = np.array([5, 2, 3])\n", "e = np.random.normal(size=nsample)\n", "y = np.dot(X, beta) + e\n", "model = sm.OLS(y, X)\n", "results = model.fit()\n", "results.params" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: y R-squared: 1.000
Model: OLS Adj. R-squared: 1.000
Method: Least Squares F-statistic: 1.875e+05
Date: Tue, 10 Nov 2020 Prob (F-statistic): 2.00e-92
Time: 21:38:35 Log-Likelihood: -75.079
No. Observations: 50 AIC: 156.2
Df Residuals: 47 BIC: 161.9
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
const 4.7206 0.457 10.332 0.000 3.801 5.640
x1 2.0984 0.211 9.931 0.000 1.673 2.524
x2 2.9893 0.020 146.287 0.000 2.948 3.030
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 1.967 Durbin-Watson: 1.919
Prob(Omnibus): 0.374 Jarque-Bera (JB): 1.894
Skew: 0.409 Prob(JB): 0.388
Kurtosis: 2.510 Cond. No. 142.


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: y R-squared: 1.000\n", "Model: OLS Adj. R-squared: 1.000\n", "Method: Least Squares F-statistic: 1.875e+05\n", "Date: Tue, 10 Nov 2020 Prob (F-statistic): 2.00e-92\n", "Time: 21:38:35 Log-Likelihood: -75.079\n", "No. Observations: 50 AIC: 156.2\n", "Df Residuals: 47 BIC: 161.9\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "const 4.7206 0.457 10.332 0.000 3.801 5.640\n", "x1 2.0984 0.211 9.931 0.000 1.673 2.524\n", "x2 2.9893 0.020 146.287 0.000 2.948 3.030\n", "==============================================================================\n", "Omnibus: 1.967 Durbin-Watson: 1.919\n", "Prob(Omnibus): 0.374 Jarque-Bera (JB): 1.894\n", "Skew: 0.409 Prob(JB): 0.388\n", "Kurtosis: 2.510 Cond. No. 142.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**R方值已经非常精确了**" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFlCAYAAAA+t0u5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de1jUZf7/8efNMCqlRh4qwWNmhkWeqDC3NG2zzJSs1LLTbr91D53bLM36dtrCos2yLcu143qujOywUWqWaZboqKhEqRkymFlJaqIwM/fvDwYWFQM5fWaG1+O6uGBuPjPzZq5tX9735z4Yay0iIiISGqKcLkBERET+R8EsIiISQhTMIiIiIUTBLCIiEkIUzCIiIiFEwSwiIhJCop0uAKBVq1a2Y8eOTpchIiJSb1auXPmjtbb1we0hEcwdO3YkMzPT6TJERETqjTHmu4raNZQtIiISQhTMIiIiIUTBLCIiEkJC4h5zRYqLi8nLy2Pfvn1Ol1KnmjRpQtu2bXG73U6XIiIiISBkgzkvL49mzZrRsWNHjDFOl1MnrLX89NNP5OXl0alTJ6fLERGREBCyQ9n79u2jZcuWERvKAMYYWrZsGfGjAiIiUnUhG8xARIdyqYbwN4qISNWFdDCHmgceeIAnnnjisL9PT09nw4YN9ViRiIhEmogJ5nSPl74TF9Fp3Hv0nbiIdI+3/mtQMIuISA1FRDCne7yMn5eFt6AQC3gLChk/L6tWwvmRRx6ha9eunH/++eTk5ADw73//mzPOOIPu3btz2WWXsXfvXpYtW8b8+fMZO3YsPXr0YNOmTRVeJyIi8lsiIpjTMnIoLPYf0FZY7CctI6dGr7ty5Upmz56Nx+Nh3rx5rFixAoDhw4ezYsUK1qxZQ0JCAi+++CJnn302Q4cOJS0tjdWrV9O5c+cKrxMRkTDz+eeQmlryvR6E7HKpI5FfUHhE7VW1ZMkSLr30Uo466igAhg4dCsC6deu49957KSgoYM+ePQwaNKjC51f1OhERCU2fvDKfs8aMINpXjC/azRdT59Lv+qF1+p4R0WOOi405ovYjUdGs6euvv55//etfZGVlcf/99x92uVNVrxMRkdCT7vGy8rW3aFS8n2gbINpXzMrX3qrzOUwREcxjB3Ulxu06oC3G7WLsoK41et1zzz2Xt956i8LCQnbv3s0777wDwO7du2nTpg3FxcXMmDGj7PpmzZqxe/fusseHu05EREJfWkYOpmg/UUAAQ7Erms/iT63xbdLKRMRQdkrPeKDkQ8wvKCQuNoaxg7qWtVdXr169GDlyJD169KBDhw6cc845ADz88MOcddZZdOjQgcTExLIwHjVqFH/605+YPHkyb7zxxmGvExGR0JdfUEiXH7eyoVVH3ks4h887nM6q+ARMDW+TVsZYa+v0DaoiKSnJHnwec3Z2NgkJCQ5VVL8a0t8qIhIu+k5cRP7OX2m2fy+7mjQta4+PjWHpuAE1fn1jzEprbdLB7RExlC0iIlKrpk/nvh7NaNLIfUAo18Zt0soomEVERMpbuhSuu44L336J1OGJxMfGYCjpKacOT6zxbdLKRMQ9ZhERkVrxyy8wejR06ACPPUZK8+Z1HsQHUzCLiIiU+tvfIC8PliyB5s0dKUHBLCIiDVa6x1u2oufq3C94eNZMePBB6NPHsZoUzCIi0iCVnrNQuqVz+vGn0arftXS86DqGOViXJn9VIi8vj2HDhtGlSxc6d+7MrbfeSlFREYsXL2bIkCGHXP/uu+/Ss2dPunfvTrdu3XjhhRccqFpERCpTes5C0tb13Lx0Fl1+zGVS8ggeX7jJ0boUzL/BWsvw4cNJSUnhm2++4euvv2bPnj1MmDChwuuLi4sZM2YM77zzDmvWrMHj8dC/f//6LVpERKokv6CQXt5sZs2+hzs+m8GM2RPo5c2u8TkLNRVZwVzLJ4AsWrSIJk2a8Ic//AEAl8vFpEmTeOmllyo8wnH37t34fD5atmwJQOPGjenatW7Xu4mISPXExcZwxdoFRAf8GMDt95Gcm1Ur5yzURPjcY66o5zliRMkMur17oW9fWLsWAgGIioLTT4dbb4Xrr4cff4TLLz/wuYsXV/qW69evp3fv3ge0NW/enPbt27Nx48ZDrm/RogVDhw6lQ4cODBw4kCFDhnDllVcSFRVZ//4REYkE95x1HGc9vBwL+E0Uxa5oPCf2qPMNRCoTPsFcmV9+KQllKPn+yy81fklrbYWnSx2uHWDatGlkZWWxYMECnnjiCT766CNeeeWVGtciIiK1yFoufvpeAkV7Sb1sLO4d29nYLYmRY4bX+7rlg4VPMP9WD/eoo2DGDBg4EIqKoFGjksel091btapSD/lgp556Km+++eYBbbt27WLr1q107tz5sM9LTEwkMTGRa665hk6dOimYRURCzc6dsGULUY9NZMIddzhdzQEiZ4y1Tx9YuBAefrjkey2sQRs4cCB79+7ltddeA8Dv9/P3v/+d66+/nqOOOuqQ6/fs2cPicv8AWL16NR06dKhxHSIiUstatIAvvoDbbnO6kkNETjBDSRiPH19rC8ONMbz11lu8/vrrdOnShZNPPpkmTZrw6KOPArBw4ULatm1b9uXxeHj88cfp2rUrPXr04P7771dvWUQklPz6K9x1V8ntzsaNS+YkhZjwGcp2SLt27XjnnXcOae/fvz+FhYdOqS89s1lERELQbbfBiy/C4MEVTyoOAaH3TwUREZG6MHcuTJsG48aFbChDFYLZGNPEGPOlMWaNMWa9MebBYHsnY8wXxphvjDFzjDGNgu2Ng483Bn/fsW7/BBERkUps2QJjxkBycsle2CGsKj3m/cAAa213oAdwoTEmGXgMmGSt7QLsBG4IXn8DsNNaexIwKXidiIhIvUv3eBnz12fY2Pt37N27nw8nPAlut9Nl/aZKg9mW2BN86A5+WWAA8Eaw/VUgJfjzsOBjgr8faA636Lfy967O08JKQ/gbRUSckO7xMmfyXJ6edicdd+bjCvh5+c3lpHu8Tpf2m6p0j9kY4zLGrAZ+AD4CNgEF1lpf8JI8oHRFdjywFSD4+1+AlhW85hhjTKYxJnPHjh2HvGeTJk346aefIjq4rLX89NNPNGnSxOlSREQiTlpGDv03fIbb7yPaWlwBPz03ryYtI8fp0n5TlWZlW2v9QA9jTCzwFpBQ0WXB7xX1jg9JV2vtVGAqQFJS0iG/b9u2LXl5eVQU2pGkSZMmtG3b1ukyREQiTlFePiPWfgiAL7jl5vL2iY4fUlGZI1ouZa0tMMYsBpKBWGNMdLBX3BbID16WB7QD8owx0cAxwM9HWpjb7aZTp05H+jQRERHw+5nywZM08RVzx8W303bXDpa3T2RVfALxDh9SUZlKg9kY0xooDoZyDHA+JRO6PgYuB2YD1wFvB58yP/j48+DvF9lIHo8WEZHQ88gjJG3yMGHIbcw/9byy5hi3y/FDKipTlXvMbYCPjTFrgRXAR9bad4G7gTuMMRspuYf8YvD6F4GWwfY7gHG1X7aIiMhhLF5csiRq9GjOePDvxMfGYID42BhShyc6fkhFZSrtMVtr1wI9K2jfDJxZQfs+4IpaqU5ERORINWoE550HU6aQ0qwZKb3Cax6PtuQUEZHIcvbZsGCB01VUm7bkFBGRyPD44yUHVPj9TldSIwpmEREJf0uXwj33lGy9GYInRh0JDWWLiEhYSvd4ScvIYe+27fz31dto1qYtR//731C9zSZDRnj/s0JERBqkdI+X8fOyOH7dKuZOv4uWu37iugvuIH3znsqfHOIUzCIiEnbSMnJI2LKOmbPv4aSf88BAoKgo5LfbrAoFs4iIhJ38gkKSc7OIDvgxgLGW5NyskN9usyoUzCIiEna6Re/HWEuxK/qAfbDjQny7zarQ5C8REQkvfj8vLZxM7Ipl3Dj0bk7+KZfl7RPJ7ngaqSG+3WZVKJhFRCS8PPooxy//BM+9j5F9dBILCwqJi40hdVDXkN9usyoUzCIiEj4WLoT774err6bnQ2NZGuZLoyqie8wiIhIedu2Cq66CU06BKVPCfr3y4ajHLCIi4aF5c3j6aTj9dGja1Olq6oyCWUREQt+2bdCmDYwa5XQldU5D2SIiEpLSPV7G/PUZ3jxtIEXt2rNk2htOl1Qv1GMWEZGQk+7xMmfyXF7+z9009hdjMUz9dDM/9fZGxMzr36Ies4iIhJy0jBySNq6isb8YAwSMIfG79RGx5WZlFMwiIhJy8gsKOW37NxjAjynb2SsSttysjIayRUQk5MQd04RNLdvxZqOj2NSyHcvbJ7IqPoH4CNhyszIKZhERCS3WMvbCUxi/9wYKi3xl65Vj3C7GRsCWm5XRULaIiISOnTuhXz9S9uWSOjyR+GOPwgDxsTGkDk+M+IlfoB6ziIg4JN3jJS0jh/zgXtdjf9+FlAf+Bp9/DtaS0jO+QQTxwRTMIiJS79I9XsbPy6Kw2A+At6CQb++8Dxa/C5Mnw9lnO1yhczSULSIi9S4tI6cslAH6blnNLZ/8hw+7D4CbbnKwMucpmEVEpN4dvOxp+PpFbGrRltsH/DViD6eoKg1li4hIvYuLjcFbUEgvbzbJuVnM6H4h3/WPI/a4Fk6X5jgFs4iI1Luxg7oyZ/JcXps5DmMDFLvc/PHqiYwc1N/p0hynYBYRkXqX0jOeM7YuwB0I3mf2+7i36Q+c2gBnYR9MwSwiIvUvM5P4Tz+CqCgwhuhGjTj1qqFOVxUSFMwiIlK/duyA4cNLzld+/nlYvRr694c+fZyuLCQomEVEpH5NmAA//ABLl0Lv3nDRRU5XFFK0XEpEROrXE0/A+++XhLIcQsEsIiL1Y9ky2LsXmjeHAQOcriZkKZhFRKTurV0L558Pt9/udCUhT8EsIiJ16+ef4dJLITYWHnzQ6WpCniZ/iYhI3fH7YfRo2LoVPvkETjjB6YpCnnrMIiJSJ9I9Xt5OHgoffMD0PpeS3qS90yWFhUqD2RjTzhjzsTEm2xiz3hhza7D9AWOM1xizOvg1uNxzxhtjNhpjcowxg+ryDxARkdCT7vEyZ/JcBnk+IgBctiydOZPnku7xOl1ayKtKj9kH/N1amwAkAzcaY7oFfzfJWtsj+PU+QPB3o4BTgQuB54wxrjqoXUREQtS0Nz6n5+bVRAf8RAFuv4+em1eTlpHjdGkhr9JgttZus9auCv68G8gGfmsz02HAbGvtfmvtt8BG4MzaKFZERMLATz/x3LM30W37Jopd0fhMFMWuaJa3TzzkuEc51BFN/jLGdAR6Al8AfYGbjDHXApmU9Kp3UhLay8s9LY/fDnIREYkUPh+MGMEJe37mljOH89IZKSTnZrG8fSKr4hOIj41xusKQV+VgNsY0Bd4EbrPW7jLGTAEeBmzw+z+BPwIVnXBtK3i9McAYgPbtNSFARCQi3HknLFpE1oOTyCk+hcJiP6viEwCIcbsYO6irwwWGvirNyjbGuCkJ5RnW2nkA1trt1lq/tTYA/Jv/DVfnAe3KPb0tkH/wa1prp1prk6y1Sa1bt67J3yAiIqHglVfg6afhttvo/X+3kTo8kfjYGAwQHxtD6vBEUnSsY6Uq7TEbYwzwIpBtrX2yXHsba+224MNLgXXBn+cDM40xTwJxQBfgy1qtWkREQk+7djBiBKSlASVnLiuIj1xVhrL7AtcAWcaY1cG2e4ArjTE9KBmm3gL8GcBau94YMxfYQMmM7huttf7aLlxERJyT7vGSlpFDfkEh7Zq5uWPwqaQMHAgDBzpdWtgz1h5y+7feJSUl2czMTKfLEBGRKkj3eBk/L4vCYj+NfUXMnHUPC7r9jq6PP6Ae8hEwxqy01iYd3K6dv0RE5IikZeRQWOynV142b0wfS+/8r9jc9DitUa4l2itbRESOSH5BIb282cyZNQ53wE9xlIsdTY/VGuVaoh6ziIgckbjYGK5c/QHRgZLpQ8ZaknOziNMa5VqhYBYRkSMydlBXCpq1wGLKdvXynNhDa5RriYayRUTkiKT0jCf92Unc+HxfOn21io3dkhg5ZrgmftUSBbOIiFRNcTFccQXccAMpl1xCygu3Ol1RRNJQtoiIVM2tt8Lbb8POnU5XEtEUzCIiUrlnn4UpU+Cuu+Daa52uJqIpmEVE5Ld99FFJb/mSS+DRR52uJuIpmEVE5LdlZEBCAsyYAS6X09VEPE3+EhGRQ6R7vLw/dR4nbchkY7ckhj47lyHNmjldVoOgYBYRkQOke7y8MWkWr0y/G4OlaNls/rhvIr5mzbUkqh5oKFtERA6Q9sFX3LlgKtE2gMta3H4fPTev1l7Y9UQ9ZhEROcCFC2bTY9s3FEe5MNZS7IpmeftE7YVdTxTMIiLyP/PnM+HjF3n/5LOZdkYKyVvXsbx9IqviE4jXXtj1QsEsIiIlrIXnn+eXhNOZcMlYduJmVdtuAMS4XdoLu54omEVEpIQx8NZbHLtrF/fnFZGWkUN+QSFxsTGMHdRVE7/qiYJZRKSh+/VXuPtueOghaNECWrcmpTUKYodoVraISEPm98Po0SXbba5c6XQ1goJZRKRBSvd4GfPXZ/iyY3d4+23W3vkA/P73TpclKJhFRBqcdI+XOZPn8q+pd3Bm3np8JorUH5qS7vE6XZqgYBYRaXDSMnJI2rgKd8BX1qYNREKHgllEpIHJLyhkacfu7I9uhM9EaQOREKNZ2SIiDUleHg9+MYMHzxjFVaMeITk3SxuIhBgFs4hIQ7FrF1x8MVdu2szMbgNZFZ/AqvgEQBuIhBINZYuINATFxXDFFbB+Pe55b/KX/3ch8bExGCA+NobU4Ylatxwi1GMWEYl01sLf/gYffgjTpsEFF5CCNhAJVeoxi4hEuo0bYeZMmDABbrjB6WqkEuoxi4hEui5dYM0a6NzZ6UqkChTMIiIRJt3j5f2p8zhv2XsEYo/h6KeeJKXnSU6XJVWkYBYRiSClu3q9/J+7aewvxmK4blJvuH207imHCd1jFhGJIGkZOZyTvYzG/mIMEDCGxO/Wa1evMKJgFhGJILu3/8jgnM8AtKtXmNJQtohIBEn5fi1xu37kwQF/4ijffu3qFYYUzCIiEaTXXX/lolZd2HR0q7I27eoVXjSULSISCe6/Hz75hJSe8dz8x/O1q1cYU49ZRCTcPfkkPPQQ7NkD/fqR0jNeQRzGFMwiImEm3eMlLSOH/IJCrvnucx6a/Qhcfjk8/rjTpUktqHQo2xjTzhjzsTEm2xiz3hhza7C9hTHmI2PMN8HvxwbbjTFmsjFmozFmrTGmV13/ESIiDUW6x8v4eVkcv34Vqf+dzH1zJrKi/WnM//tEcLmcLk9qQVXuMfuAv1trE4Bk4EZjTDdgHLDQWtsFWBh8DHAR0CX4NQaYUutVi4g0UGkZOSRsWceM2RMYsfZDXDbA08kjeGzxd06XJrWk0mC21m6z1q4K/rwbyAbigWHAq8HLXgVSgj8PA16zJZYDscaYNrVeuYhIA5RfUEhybhZuv48oSjYQOf37jVqnHEGOaFa2MaYj0BP4AjjeWrsNSsIbOC54WTywtdzT8oJtB7/WGGNMpjEmc8eOHUdeuYhIA3Sqq5DzNq3AF+U6YAOROK1TjhhVnvxljGkKvAncZq3dZYw57KUVtNlDGqydCkwFSEpKOuT3IiJykF27mD7vQRr9sJkHzv8zLQp3sbx9ItkdTyNV65QjRpWC2RjjpiSUZ1hr5wWbtxtj2lhrtwWHqn8ItucB7co9vS2QX1sFi4g0SPv2QUoKsRu/Ytmkl1iyJ578gkLiYmNIHdRVy6MiSKXBbEq6xi8C2dbaJ8v9aj5wHTAx+P3tcu03GWNmA2cBv5QOeYuISOXKL4eKi41h7PknkZJ6O3z8MUyfztmjR7PU6SKlzlSlx9wXuAbIMsasDrbdQ0kgzzXG3ADkAlcEf/c+MBjYCOwF/lCrFYuIRLDS5VCFxX4AvAWFpM7+gv4bviH26adh9GiHK5S6Zqx1/vZuUlKSzczMdLoMERHH9Z24CG+5GdbGBrAmio5NXSy+90IHK5PaZoxZaa1NOrhdO3+JiISQ0mVPvbzZ3LhsLrH7dnHVqEf5bk9jhyuT+qJgFhEJIXGxMRy/fhWzZ46nUcCH30Rx2vcb+T7xkI6VRCidLiUiEkLGDurKtWv+izvgA0rWmvbdlq1jGxsQ9ZhFREJISsHX+LOXYI3Bj8EX7ab3tZfST8uhGgwFs4hIKDnmGFzn/A7uuosoj4fo/v3p16eP01VJPVIwi4iEgh9/hFatoFcvWLgQjIELNQu7IdI9ZhERp23aBImJ/ztP+fBbHksDoB6ziEg9K7+z1+lmDzNfG8vRxcUwZIjTpUkIUDCLiNSj8jt7tdj7C0/MHIfd/SOLX3qD/t26OV2ehAAFs4hIPUrLyKGw2E/vret55p3Habn3F64Z+Q+8W5to/2sBFMwiIvUqv6CQXt5sps+9j0a+YnwuF74oV9mOXyKa/CUiUo86Hu3i0nWLcPt9uLC4AgGSc7OIi41xujQJEeoxi4jUl6IiZn30T1qs/RifKxr8Popd0XhO7KGdvaSMgllEpD74fDB6NCcsWcDq8Y/y3M6mnLQhk43dkhg5Zjgp2tlLghTMIiJ1ze+HP/wB3ngDnnySHrffzlSna5KQpXvMIiJ1bdYsmD4dHnkEbr/d6WokxKnHLCJS10aPhhYtYPBgpyuRMKAes4hIXbAWHnsMNm4s2WJToSxVpGAWEalF6R4vY/76DEs79oBx48iZ+IzTJUmYUTCLiNSSdI+XOZPn8uzU2+mbuxafieIBX3vSPV6nS5MwomAWEaklaRk5/OWT6bgD/rK2nt+uJS0jx8GqJNxo8peISC3Z8dMujt/9Ez5T0ucpdkWzvH2ittuUI6JgFhGpDT4frVs2Z/g1T9Bt+2bOzFvP8vaJrIpPIF7bbcoR0FC2iEhNTZoEF1zAuHPbY49uSma7U3muzwhWxScQ43Zpu005IgpmEZGaePppuOMOaNGCS87oQOrwROJjYzBAfGwMqcMTtd2mHBENZYuIHIF0j5e0jBzyCwq5eUMGd7zzDFx6acnuXm43KT3jFcRSIwpmEZEqSvd4GT8vi4Qt67j3i3lc9M3nLOjah1/H/ZNhbrfT5UmEUDCLiFRRWkYOCVvWMWP2BBr5i/GZKF7oPZT8Rd8y7MxOTpcnEUL3mEVEqii/oJDLshbi9vtwWQvAGXnZWg4ltUrBLCJSRbeve4/Raz4gYKLwmaiydcpxWg4ltUhD2SIiVfH449zy3hQyTunLtF6XcEbeBpa3TyS742mkajmU1CIFs4hIZR55BO69F0aOZN8dE8lftJkp7U4jLjaG1EFdNQtbapWCWUTkt2Rnw/33l5yp/MorDIuOZtiZHZ2uSiKYgllE5LckJMBnn8EZZ4DL5XQ10gBo8peIyMGWLYP+/WHixJLHyckKZak36jGLiJTzyStv0/eGy4gO+PEt+YylJ3Sj3/VDnS5LGhD1mEVEgtJXbiV2wl1El56nbC0rX3uLdI/X2cKkQak0mI0xLxljfjDGrCvX9oAxxmuMWR38Glzud+ONMRuNMTnGmEF1VbiISK3y+2n8pxvonv81xVGusnXKn8WfSlpGjtPVSQNSlaHsV4B/Aa8d1D7JWvtE+QZjTDdgFHAqEAcsMMacbK3110KtIiJ1JyoKr+sonjjnapa1707y1qyy85SNdvaSelRpMFtrPzXGdKzi6w0DZltr9wPfGmM2AmcCn1e7QhGRulRYCPn50LkzLw+/Ge8v+wBY1Tah7BLt7CX1qSb3mG8yxqwNDnUfG2yLB7aWuyYv2HYIY8wYY0ymMSZzx44dNShDRKSadu+GwYOhXz/49VfGXngKMe4DZ1/HuF2M1c5eUo+qG8xTgM5AD2Ab8M9gu6ngWlvRC1hrp1prk6y1Sa1bt65mGSIi1bRzJ1xwASxZAo89BkcfTUrPeFKHJxIfG4MB4mNjSB2eqJ29pF5Va7mUtXZ76c/GmH8D7wYf5gHtyl3aFsivdnUiIrUs3eNlydOvcdcbabQo3MWqx5/nrNGjy36f0jNeQSyOqlaP2RjTptzDS4HSGdvzgVHGmMbGmE5AF+DLmpUoIlI70j1e5kyey2Ov3ctxv+4kYAxPrdut5VASUqqyXGoWJZO3uhpj8owxNwCPG2OyjDFrgfOA2wGsteuBucAG4APgRs3IFpFQkZaRQ8/Nq4GS+26uQICem1drOZSElKrMyr6yguYXf+P6R4BHalKUiEitW72aB168h5d6X0KxKxr8vrLzlPO1HEpCiLbkFJHIt2QJDBlComlCXuwJjB71CMm5/1unHK/lUBJCFMwiEtneew8uvxw6dGDNk6/x4+c/szX2BFbFl6xT1nIoCTUKZhGJGOkeL2kZOeQXFBIXG8OT0Zs4a/zfoEcPeP99BrVuTWqbA68ZO6irZmFLSDHWVrjMuF4lJSXZzMxMp8sQkTCW7vEyfl4WhcX/m2968u7tvLRpPm1f/w80a+ZgdSKHMsastNYmHdyu06VEJCKkZeRQWOynV142T777T3rlZfN1s+MZ2f8WhbKEFQ1li0hEyC8opPfW9cyZNZ5oG2BI9hJGXZWKh4TKnywSQtRjFpGI0PFoF49++CzRNgBAlA2QnJulAygk7KjHLCLhr6CA19MfoNWPuRRHuTDWUuyKxnNiD824lrCjYBaR8JeVRaucdax49F/8OzfASRsy2dgtiZFjhmvGtYQdBbOIhK9ffoFjjoFzzoEtWzijVSvOcLomkRrSPWYRCRvpHi99Jy6i07j3+PNfn2F/x07wxhslv2zVytniRGqJeswiEhZK1yknbFnHHZ73GfLVErYdczxfN47jAqeLE6lFCmYRCQtpGTkkbFnH7JnjcQd8WAwPDPgT36zfxwWXOF2dSO3RULaIhIX8gkIuXbeIRgEfBggYQ7cfNutkKIk4CmYRCQtxsTG8ddoAiqKi8ZmosiMbtU5ZIo2GskUktP38M1xzDf+4+hb+9utpjLoqtezIxuyOp5GqdcoSYdRjFpHQtXkznH02LFjAeY1/JXV4IttP7cWUPiPYfmovUocnap2yRBz1mEUkNH3xBVxyCfh8sGABnHMOKaAgloinYBaRkJHu8fL+1Hn0/TKDUWs+xNcmjgZHt6UAABoASURBVKOXZEBXDVdLw6FgFpGQkO7xMmfyXF6aPg6334c1htv7/pHBe5uS4nRxIvVI95hFJCQ8+f4Gblr0Cm6/j2gbwFjLSds2kZaR43RpIvVKwSwiziso4B/T7qZv7loCJuqA5VBapywNjYayRcRZmzbBkCGcnfsNYy+6hU0t25Uth1oVn0C81ilLA6NgFpF6ke7xkpaRQ35BIXGxMYwd1JWUwPcwaBBYy/LnZ/FubjMKi/2sik8AIMbt0nnK0uBoKFtE6lzpARTegkIs4C0oZPy8LN7b1RiSk+GLL/jdn64gdXgi8bExGCA+NkbrlKVBUo9ZROpcWkYOhcV+enmzSc5dS7P9e5n0u6t59IsfuPjdd8uuS+kZryCWBk/BLCJ1Lr+gkF7ebGbMnkBjXxFRwH6Xm6fPudrp0kRCjoayRaTOxcXGcMHXn9MkGMp+DEXRjXQAhUgF1GMWkTqXevwueqz5ACgNZTeeE3toYpdIBRTMIlLnzj3rZH7ucjL/d9oQmuXnsrFbEiPHDNf9ZJEKKJhFpFYcvBzqrgGdGJbzGVx9NXTrRou1q3jYGKfLFAl5CmYRqbHS5VCFxX4A9nvzaXvZzZC3ATp3Ljm6UaEsUiUKZhGpsfLLoYZuWMzFXy2laVEh/3flvTx09tlOlycSVhTMIlJjpcuhZs8cjzvgwwJjL7qVee2Tecjp4kTCjJZLiUiNxcXGkJybhSvgxwABE8Xxv+7UciiRalAwi0jNfP89k6K+wXNiD4qi3WUnQ2k5lEj1aChbRKrv88/h8ss5c88ern57Kbc2ieakDZlaDiVSA5UGszHmJWAI8IO19rRgWwtgDtAR2AKMsNbuNMYY4GlgMLAXuN5au6puShcRx1gLU6bAbbdB+/bw3/9y8emncXH/05yuTCTsVWUo+xXgwoPaxgELrbVdgIXBxwAXAV2CX2OAKbVTpog4Ld3jpe/ERXS6+13eOXMw3HgjXHABrFgBp5/udHkiEaPSYLbWfgr8fFDzMODV4M+vAinl2l+zJZYDscaYNrVVrIg4o3Sd8vHrV/HX5a9T4DM8c+7VpD84BY491unyRCJKde8xH2+t3QZgrd1mjDku2B4PbC13XV6wbVv1SxQRp6Vl5HD1Z3O5c8l/cAUCFLuiGT3qEbZ/9A0pvds5XZ5IRKntWdkVbe1jK7zQmDHGmExjTOaOHTtquQwRqTU+H9ekP8eExS/TyO8j2gZw+30k52aRX1DodHUiEae6wby9dIg6+P2HYHseUP6fz22B/IpewFo71VqbZK1Nat26dTXLEJE65fXCgAH85Ys3yehyFvuiG5Uth1rePlHrlEXqQHWHsucD1wETg9/fLtd+kzFmNnAW8EvpkLeIhLaDD6H4v+7NGHTdxbB3L5n/mMxt+7qQsGUdyblZLG+fSHbH00jVOmWRWmesrXCk+X8XGDML6A+0ArYD9wPpwFygPZALXGGt/Tm4XOpflMzi3gv8wVqbWVkRSUlJNjOz0stEpI4cfAgFQEx0FG/nv8fJY2+EhIRDgnvsoK5apyxSA8aYldbapEPaKwvm+qBgFnFW34mL8BYUct7GLxn/ySs8dfYo3k84l/jYGJaOG+B0eSIR6XDBrJ2/RIT8gkKuXfkuDy54HoCn3nuS75u3xkOCw5WJNDzaK1ukofP5uP/L2TwQDGUDRAUCJOdmaXKXiAMUzCIN3aRJXP/xdD456YwDZl3rEAoRZ2goW6QBqHDiVudm0Lw53HQTdOnCLx3O4Nap83QIhYjDNPlLJMIdPOO6cfF+Hlz8IkN+zKbpujXQrJnDFYo0TJr8JdJApWXkUFjsp5c3m4uzlzBw45d0/OV7Zpw7gtGNGztdnogcRMEsEuHyCwrplZfN7FnjcQd8ADw04P/x8hkpjG7UyOHqRORgmvwlEuHiYmPok7sWd8CHAfwmiia+Is24FglRCmaRCFB2VvK49+g7cRHpHm/JL957jwlntGRV556acS0SJjSULRLmDp7c5S0o5KHZK+j+yDg6vTmdwbffTtEtf+fWJtGacS0SBhTMImGudHJXqdO3fc1T7zxBh4JtcNdd8NBDpDRuTMqUmx2sUkSqSsEsEuZKz0Tu5c3mmlXvcsmGT9nerBVXjXqE2Y+Nd7g6ETlSCmaRMBcXG8Px61YxY84E3H4fGMPdF93M1tPPcro0EakGBbNIOAsEeP7HJTR/73ncfh/RNoDPRNFrx2YuH/RHp6sTkWrQrGyRcLV5M5x3Hon/fICm7dvic0XjM1H4ot30vvZSTe4SCVPqMYuEm0AAXngBxo4FlwtefpmW110Hy5fD4sVE9+9Pvz59nK5SRKpJwSwSBtI9Xt4PHjDx3cndmfjRszTr2xemTYN27Uou6tOn5EtEwpqCWSTEpXu8zHl6Dq/+526ibIDiZbO5ffg9XHj3DaS0a+t0eSJSy3SPWSTEzZj1MY++/QSNAiWTu9x+Hyd+v5m0D792ujQRqQPqMYuEKr8fJk/m1afGg7UURUWX9Jhd0Sxvn1i2fllEIouCWSRU3X47PPMMnq5nceeAv9Bm948k52axvH0iq+ITiNchFCIRScEsEkqKiuDXX+HYY+HmmyE5mR2nnEvBW+vY1rw1q+ITAIhxu3QIhUiEUjCLOKx0xvU5yz9gwHerMD16ELfov9ClC3TpQgqAMaRl5JBfUEhcbAxjB3XVOmWRCKVgFnFQusfL20+8yr9n3YfLBrBAWlRbunq8BwRvSs94BbFIA6FgFqlj6R7vYXu77057m8lv/AOXDWAAv4nCFBeRlpGjIBZpoLRcSqQOlZ6V7C0oxFJyVvL4eVmkr9wKwCpXLBtbtKPI5cZnojTjWkTUYxapSweflez2F/OHz+fSYcZq+GYVMXEncOl1T9LLm60Z1yICKJhF6lT5s5Ivy1rAOd+upv2u7fz35LNhzx7GDurK+HlZrIpP0IxrEQEUzCJ1Ki42hk6eZbzy+v1E2wAB4B/n/ZH/XjCai445hpSexwBoxrWIlFEwi9ShsYO68t0Hr+CyAQACJoqjCBzQI9aMaxEpT5O/ROrCokUwaBApXY6hx/WXs9/dWGcli0iVqMcsUkPll0P1ZBfPrJxB/IL34MQTYcsW+l0/FLp+rLOSRaRKFMwiNVC6HOrUzWtJXTaLs7auJ2Ci2PC3sXT750PQpEnJhTorWUSqSMEsUgNpGTkkbFnHf+beRxNfEQETxV9TxrK+3UCWloayiMgRUDCLVNeXXzLxhTtZd/yJuP0+DGCBk37K4yNtECIi1aRgFqnEwVtq/l/3Zgya8TTMmMGpTWP56KSzKHZFg99XtnNXnDYIEZFqUjCL/IbSe8gJW9YxLDeLDju30e++T/BHGVz33MOyi67l9YxvWXdC57Kdu7I7nkaqNggRkWpSMIv8hrSMHBK+XceMORNw+30ALGt/Ok+NvIt5j1zJEMB3dFPSMhoxJT6BuNgYUrVBiIjUQI2C2RizBdgN+AGftTbJGNMCmAN0BLYAI6y1O2tWpogDrOXUFR/zyAfP0MhXjAuLz0SxvMPpeEzzssu0QYiI1Kba2GDkPGttD2ttUvDxOGChtbYLsDD4WCS8fPYZ/O53TJ33D/a73Phc0Qec/qR7yCJSV+piKHsY0D/486vAYuDuOngfkboxejTMnAlt2uCZMJFrA6fSZetXuocsIvXCWGur/2RjvgV2UrJK5AVr7VRjTIG1NrbcNTuttcdW8NwxwBiA9u3b9/7uu++qXYdIdaV7vLw/dR49PZ/i6Xkug8cMJ2XxXNi/H265BY466pBZ2TpkQkRqgzFmZbnR5jI17TH3tdbmG2OOAz4yxnxV1Sdaa6cCUwGSkpKq/68DkWpK93jJmDiN515/CJcNULwinev2+eCWEQcEr+4hi0h9qtE9ZmttfvD7D8BbwJnAdmNMG4Dg9x9qWqRIrdu+ncJbbmPyGw/jsgEMEGUD9Ny8mrSMHKerE5EGrNrBbIw52hjTrPRn4AJgHTAfuC542XXA2zUtUqRWWQsDBzJi6Tw+7dizZHJXuYld+dq1S0QcVJOh7OOBt4wxpa8z01r7gTFmBTDXGHMDkAtcUfMyRaqn9B7yaWuW0qF4N2by0wztcxI8+yyj529hubsVvbzZZRO7VsUnEK8Z1yLioGoHs7V2M9C9gvafgIE1KUqkNqR7vLz3+Es8N/dBogN+AB6/L45A2gOk9OvHqOYnsWZeFqviE1gVnwBAjNvFWM24FhEH1cY6ZpHQU1RE4S238ezcB3EH/BjAb6IwxUVl95BTesaTOjyR+NgYDBAfG0Pq8ERN9BIRR2lLToksv/wCxxwDbjedt2SztEN3+uRmER3wV3gPWTOuRSTUKJglrJXeQz4rcyHdduWT5M3GvXkTHHccd/5lErm7i3QPWUTCioJZwla6x0vmfWlMef8ZomwAgPe7nYtZt53BA47jjsHdGK97yCISZhTMErIq23HrtdmfMu+9p7GAAXwmivWtOvL2lzsYPICya7Vrl4iEEwWzhKTSc5ALi0tmU3sLCnlwzgo6vv4aPfb/CP/8Jx7TnNR+13Pb0pm4/T7dQxaRiKBglpCUlpFDYbGfXt5sBmz8kta/7uT3G7+kReEuSE6GoiLiYmN4IflyVrQ7VfeQRSRiKJglJOUXFNLLm83sWeNx+30ArIjvxhP9rmXu9LvAGMYO6qp7yCIScRTMElqKiuDNN7n45220z80iKhAou3+8uHMS3sQzoGS3Od1DFpGIpGAWx5QudTppQyY7OpzMDTE/ccr8WfD999x90aXcdeI5FC+bDcH7x54TexzSG9Y9ZBGJNApmcUS6x8ucyXN5afo4GvuKMJTMrP7+dwM44aWXaDdoECPXbOPWJtGctCGTjd2SGDlmuEJYRCKeglnqzGGXO/34I9/d9yhnFBTg9vuIAgLAv8+4lFcvvomlFw0Agr3hKTc7+jeIiNQ37ZUtdaJ0uZO3oBALeHfu5a0np5N3YQrEx3Pre8+xL7oxxa5ofCaK/dGN+KDr2TpyUUQaPPWYpU6UX+503qYVXLZuIXG7f2JXk6bwl79wtTmdz5qcQGa7blrqJCJSjoJZqu2wQ9X79tFj+Uf80buBq1Z/ULbc6amzR/F88hV89c/hXO7xslJLnUREDqFglmo5ZGeunXuZ8/Qcuu/x0GnhuzxbUMBPMc1x+31E2wA+E0VRdCNato4FtNRJRORwFMxSLeWHqpNzs2hb8D1Xrf2QQncTGHUFS88ezLOf5/HizAmHXe6kpU4iIodSMMuR83oZ9NFsRq7NoOPObbgCfnxRLib3GcnUsy5j3aQr6AvsOMvLrUc11nInEZEjoGCWCpXf/GNjtySGXH8xQ1d+ALNnw2ef8X/Wsv3oY4n2+3BhIQD73I055viWZa+h5U4iIkdOwSyHKN384+Xp43D7iilaNpsb9u7ngvcep8lxreCBB1iQ2J8X56/kpenjyk52qmhnLhEROTIK5gbqsDOqN28m974nePyzeWU7crn9PnpsyeLyPz/Hu/+4DIzhfGBPx87amUtEpJYpmBug0hnVCVvWMSw3i+XtEhn/axEJzz5O1xcncwuwJfYEfFEujLVl5xyv98eUHSABGqoWEakLCuYGaPI7axi99E3GLX4Flw2w3+XmqisfZdLRnXn+ySe5YltrVkQdWzbjWpt/iIjUHwVzhDp48tbgMcNJiS2Cv/yF/y78mMb+YiwlB0dEB/wk52Yxpc8IuP1iRnu8rNPmHyIijlAwR6B0j5f5T7zKlDn34/b78H02k+v2+XD9v4u55PvveSt5KF81bsHdn75aNnFreftE4oI9Ym3+ISLiHAVzpHn4YbpMncm0vByisABEB3z03LyaiUu6c8maNTTxeJkzL4u1cSeXDVVndzyNVG3+ISLiOAVzGEr3ePnouTn87ssPadooisS2x9DxrVklv1yxgr02itcTB5Ky4RNcAX9Zj7j05Kb/9YgbMSU+gbjYGFLVIxYRCQnGWut0DSQlJdnMzEynywgtn38OixdD//7Qp09Zc+Yjz3DMU09w0o+5GMACG1u1Z8N7ixl2Ziewlr6PfYy3oLDCyVtLxw1w6A8SEZHyjDErrbVJB7erx1wHDrtGuIo+ffFNzvrLlUT7fGAMRS1bE5P5BXTowIJVWxhdtK9s4pbfRPFWt/68vejbkmA2hrGDujJek7dERMKSgrmWHXLqUkEh4+dlAf8bQi4/Y/rbU3oy4rxunNf/dDjhBJZOmcXZN15NtA0AYK0lzzThu5VbOL9DB17oMoAvj2rDjNkTDpi4VTpMXf59NHlLRCT8KJhrWempS+UVFvtJy8ghpWc8736Ww8933sNzX6bjsgH49DXMVFh710Oc/th9pOVGcf0p5zA4ZylRNkCxK5q7B93E9q/9nA/ExcawigRGj3rkN9cYa/KWiEh4anDBXNNh5sqen19QSLP9v3JhzlLO25RJsSuaYwt3s6xjdxg3gKcWbmTBF/PKhqIDwLunnMO0qJN5G1hjm3Lb0LG85h1yQPCaYI9Yw9QiIpEtooK5stCs6jDz4V6j/FaWl3/rYesxx7N4xUJObuOnW2In+POfiYuN4b/3j6B50V6gZHLW5mPjiTn6KAA2FcINw+/jX/MfKxuKfjlpKGttU6CkR+wtKDwgeEvby9epYWoRkcgUMcFcldAtHWY+eLZy6TBzWfB+m8WVm1bw3bFtWJTZjM4doklsewxp+7uTsGUdr8+4C9fBs9kHDYI//5mxg7qy9LXeXJD9GS4sfhPF2z1+z4kP3wOUBOzCLmcddii6tEdcfjj84B6xhqlFRCJXxATzb4Zut1ZQUEB+cAnRzFn34Pb7CERF8W7X3xFwuWD1VNJ6jCkJ3pnjSu7/lnfSSeRf9hTDcrMI7tuBH8N/el1MWr/rWP/k5UBJaH5y150Uj/kS6yvGF+2m97WX0i8YpJUNRatHLCLSsEVMMJeG7ozZE2jkK8YABU2acpRvP4wvAqD9wx+QnJtFI38xUYAr4OeSr5bwY/NWcHQn8jvsDQZvSfL6MbzW+2Ke+t3VrJk0grjHPmZ5+0SKot1lw9Dzu/Uj9rgWB9TS7/qh0PVjWLyY6P796VduHXJVglc9YhGRhitigjkuNobkz7Nw+324sASA3NgTWH9yT0YP6g6xsfy9dxdmbepB0bJZRPtLdsT64+hURt46kpSe8cRNXHRI8L6T0I+mJ7T+3/rgX4sOGIY+eCvLMn36HLAxSHkKXhEROZw6C2ZjzIXA04ALmGatnVhX7wUlQ8RzcnpQvGw2BEP1sQv/yshbRkAwBIcCgUYjuKVJdNmpSyPHDC8LycqCV1tZiohIXauTLTmNMS7ga+D3QB6wArjSWruhoutra0vOCo86PMLQrOlyKhERkao43JacdRXMfYAHrLWDgo/HA1hrUyu6Xntli4hIQ3O4YI6qo/eLB7aWe5wXbCtf0BhjTKYxJnPHjh11VIaIiEh4qatgNhW0HdA1t9ZOtdYmWWuTWrduXUdliIiIhJe6CuY8oF25x22B/Dp6LxERkYhRV8G8AuhijOlkjGkEjALm19F7iYiIRIw6WS5lrfUZY24CMihZLvWStXZ9XbyXiIhIJKmzdczW2veB9+vq9UVERCJRXQ1li4iISDUomEVEREKIgllERCSEKJhFRERCiIJZREQkhNTJXtlHXIQxO4DvavElWwE/1uLrNVT6HGtOn2HN6TOsOX2GNVcXn2EHa+0hW1+GRDDXNmNMZkUbg8uR0edYc/oMa06fYc3pM6y5+vwMNZQtIiISQhTMIiIiISRSg3mq0wVECH2ONafPsOb0GdacPsOaq7fPMCLvMYuIiISrSO0xi4iIhKWIC2ZjzIXGmBxjzEZjzDin6wk3xph2xpiPjTHZxpj1xphbna4pXBljXMYYjzHmXadrCVfGmFhjzBvGmK+C/5vs43RN4cYYc3vwv+V1xphZxpgmTtcU6owxLxljfjDGrCvX1sIY85Ex5pvg92Pr6v0jKpiNMS7gWeAioBtwpTGmm7NVhR0f8HdrbQKQDNyoz7DabgWynS4izD0NfGCtPQXojj7PI2KMiQduAZKstadRcgzvKGerCguvABce1DYOWGit7QIsDD6uExEVzMCZwEZr7WZrbREwGxjmcE1hxVq7zVq7Kvjzbkr+jzDe2arCjzGmLXAxMM3pWsKVMaY5cC7wIoC1tshaW+BsVWEpGogxxkQDRwH5DtcT8qy1nwI/H9Q8DHg1+POrQEpdvX+kBXM8sLXc4zwUKtVmjOkI9AS+cLaSsPQUcBcQcLqQMHYisAN4OXhLYJox5miniwon1lov8ASQC2wDfrHWfuhsVWHreGvtNijpwADH1dUbRVowmwraNO28GowxTYE3gdustbucriecGGOGAD9Ya1c6XUuYiwZ6AVOstT2BX6nD4cNIFLwPOgzoBMQBRxtjrna2KqlMpAVzHtCu3OO2aNjmiBlj3JSE8gxr7Tyn6wlDfYGhxpgtlNxOGWCMme5sSWEpD8iz1paO2LxBSVBL1Z0PfGut3WGtLQbmAWc7XFO42m6MaQMQ/P5DXb1RpAXzCqCLMaaTMaYRJZMc5jtcU1gxxhhK7ullW2ufdLqecGStHW+tbWut7UjJ/wYXWWvVSzlC1trvga3GmK7BpoHABgdLCke5QLIx5qjgf9sD0QS66poPXBf8+Trg7bp6o+i6emEnWGt9xpibgAxKZh++ZK1d73BZ4aYvcA2QZYxZHWy7x1r7voM1ScN1MzAj+A/tzcAfHK4nrFhrvzDGvAGsomTFhQftAlYpY8wsoD/QyhiTB9wPTATmGmNuoOQfPFfU2ftr5y8REZHQEWlD2SIiImFNwSwiIhJCFMwiIiIhRMEsIiISQhTMIiIiIUTBLCIiEkIUzCIiIiFEwSwiIhJC/j+XF34oVq7VcwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "y_fitted = results.fittedvalues\n", "fig, ax = plt.subplots(figsize=(8,6))\n", "ax.plot(x,y,'o',label='data') # 原始数据\n", "ax.plot(x,y_fitted,'r--.',label='OLS') #拟合数据\n", "ax.legend(loc='best')\n", "plt.show() # 可以看到图已经非常精确了" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**分类变量：**\n", "假设分类变量有3个取值（a,b,c），比如考试成绩有3个等级。a就是（1,0,0），b（0,1,0），c（0,0,1），这个时候就需要3个系数β0，β1，β2，也就是β0x0+β1x1+β2x2" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n", " 0, 0, 0, 0, 0, 0])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nsample = 50\n", "groups = np.zeros(nsample, int)\n", "groups" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [1., 0., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 1., 0.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.],\n", " [0., 0., 1.]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#亚变量，转化成类onehot的形式，变成0，1的形式\n", "groups[20:40] = 1\n", "groups[40:] = 2\n", "dummy = sm.categorical(groups, drop=True)\n", "dummy" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: y R-squared: 0.994
Model: OLS Adj. R-squared: 0.994
Method: Least Squares F-statistic: 2589.
Date: Tue, 10 Nov 2020 Prob (F-statistic): 2.80e-51
Time: 21:53:52 Log-Likelihood: -74.545
No. Observations: 50 AIC: 157.1
Df Residuals: 46 BIC: 164.7
Df Model: 3
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
const 7.9672 0.635 12.551 0.000 6.689 9.245
x1 2.0457 0.073 28.010 0.000 1.899 2.193
x2 -0.0252 0.403 -0.063 0.950 -0.835 0.785
x3 2.4251 0.336 7.209 0.000 1.748 3.102
x4 5.5673 0.758 7.346 0.000 4.042 7.093
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 0.042 Durbin-Watson: 1.779
Prob(Omnibus): 0.979 Jarque-Bera (JB): 0.224
Skew: -0.028 Prob(JB): 0.894
Kurtosis: 2.677 Cond. No. 1.28e+17


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.16e-31. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: y R-squared: 0.994\n", "Model: OLS Adj. R-squared: 0.994\n", "Method: Least Squares F-statistic: 2589.\n", "Date: Tue, 10 Nov 2020 Prob (F-statistic): 2.80e-51\n", "Time: 21:53:52 Log-Likelihood: -74.545\n", "No. Observations: 50 AIC: 157.1\n", "Df Residuals: 46 BIC: 164.7\n", "Df Model: 3 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "const 7.9672 0.635 12.551 0.000 6.689 9.245\n", "x1 2.0457 0.073 28.010 0.000 1.899 2.193\n", "x2 -0.0252 0.403 -0.063 0.950 -0.835 0.785\n", "x3 2.4251 0.336 7.209 0.000 1.748 3.102\n", "x4 5.5673 0.758 7.346 0.000 4.042 7.093\n", "==============================================================================\n", "Omnibus: 0.042 Durbin-Watson: 1.779\n", "Prob(Omnibus): 0.979 Jarque-Bera (JB): 0.224\n", "Skew: -0.028 Prob(JB): 0.894\n", "Kurtosis: 2.677 Cond. No. 1.28e+17\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "[2] The smallest eigenvalue is 4.16e-31. This might indicate that there are\n", "strong multicollinearity problems or that the design matrix is singular.\n", "\"\"\"" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Y = 5+2X+2Z1+6*Z2+9*Z3，Z1、Z2、Z3分别表示取a、b、c\n", "x = np.linspace(0,20,nsample)\n", "X = np.column_stack((x, dummy))\n", "X = sm.add_constant(X)\n", "beta = [5,2,3,6,9] #即上面的假设\n", "e = np.random.normal(size=nsample)\n", "y = np.dot(X, beta) + e\n", "result = sm.OLS(y,X).fit()\n", "result.summary()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import pandas as pd\n", "from plotly.offline import init_notebook_mode,iplot # 可以交互的包，鼠标对着有结果\n", "import plotly.graph_objs as go\n", "\n", "init_notebook_mode(connected=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }