diff --git a/doc/pub/week38/ipynb/week38.ipynb b/doc/pub/week38/ipynb/week38.ipynb index 919581e58..4ab7262cd 100644 --- a/doc/pub/week38/ipynb/week38.ipynb +++ b/doc/pub/week38/ipynb/week38.ipynb @@ -238,10 +238,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", @@ -402,10 +399,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "#Model training, we compute the mean value of y and X\n", @@ -513,10 +507,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "X = X - np.mean(X,axis=0)" @@ -724,10 +715,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -891,10 +879,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -978,12 +963,138 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Beta values for own Ridge implementation\n", + "[ 3.43579948e-02 -5.43330971e-01 -3.10141414e-03 2.47116868e-01\n", + " 2.18613217e-01 1.02054837e-01 -4.25617662e-04 -5.90475506e-02\n", + " -7.68534263e-02 -6.68929213e-02 -4.24906604e-02 -1.40927184e-02\n", + " 1.11482289e-02 2.88529063e-02 3.67047975e-02 3.38135733e-02\n", + " 2.02198702e-02 -3.46383924e-03 -3.63025821e-02]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[ 3.43579948e-02 -5.43330971e-01 -3.10141413e-03 2.47116868e-01\n", + " 2.18613217e-01 1.02054837e-01 -4.25617658e-04 -5.90475506e-02\n", + " -7.68534263e-02 -6.68929213e-02 -4.24906604e-02 -1.40927184e-02\n", + " 1.11482289e-02 2.88529063e-02 3.67047975e-02 3.38135733e-02\n", + " 2.02198702e-02 -3.46383925e-03 -3.63025821e-02]\n", + "Intercept from own implementation:\n", + "1.0330308045181225\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "1.033030804518383\n", + "MSE values for own Ridge implementation\n", + "3.139255958275475e-06\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "3.139255958572018e-06\n", + "Beta values for own Ridge implementation\n", + "[-0.05807125 -0.29822833 -0.08551306 0.08156108 0.13679863 0.12333649\n", + " 0.08251519 0.03815288 0.00111756 -0.02498832 -0.04010697 -0.04566964\n", + " -0.04355837 -0.03562355 -0.02348765 -0.00848904 0.00831018 0.0260906\n", + " 0.04423486]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[-0.05807125 -0.29822833 -0.08551306 0.08156108 0.13679863 0.12333649\n", + " 0.08251519 0.03815288 0.00111756 -0.02498832 -0.04010697 -0.04566964\n", + " -0.04355837 -0.03562355 -0.02348765 -0.00848904 0.00831018 0.0260906\n", + " 0.04423486]\n", + "Intercept from own implementation:\n", + "1.0411487294305548\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "1.0411487294305266\n", + "MSE values for own Ridge implementation\n", + "1.9601304850163794e-05\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "1.9601304850085328e-05\n", + "Beta values for own Ridge implementation\n", + "[-0.1416398 -0.14021063 -0.05383795 0.01367553 0.04784395 0.05796251\n", + " 0.05447415 0.044613 0.03267527 0.02098261 0.01066519 0.00217499\n", + " -0.00440346 -0.00917248 -0.01231917 -0.01405935 -0.0146081 -0.01416528\n", + " -0.01290947]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[-0.1416398 -0.14021063 -0.05383795 0.01367553 0.04784395 0.05796251\n", + " 0.05447415 0.044613 0.03267527 0.02098261 0.01066519 0.00217499\n", + " -0.00440346 -0.00917248 -0.01231917 -0.01405935 -0.0146081 -0.01416528\n", + " -0.01290947]\n", + "Intercept from own implementation:\n", + "1.0495569966278282\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "1.0495569966278269\n", + "MSE values for own Ridge implementation\n", + "5.4959161509370406e-05\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "5.4959161509366834e-05\n", + "Beta values for own Ridge implementation\n", + "[-0.13535942 -0.08593216 -0.03568439 -0.0036367 0.01397146 0.02229529\n", + " 0.02503753 0.0245528 0.02228115 0.01908936 0.01549377 0.01179792\n", + " 0.00817631 0.00472512 0.00149311 -0.00149956 -0.00424967 -0.00676387\n", + " -0.00905423]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[-0.13535942 -0.08593216 -0.03568439 -0.0036367 0.01397146 0.02229529\n", + " 0.02503753 0.0245528 0.02228115 0.01908936 0.01549377 0.01179792\n", + " 0.00817631 0.00472512 0.00149311 -0.00149956 -0.00424967 -0.00676387\n", + " -0.00905423]\n", + "Intercept from own implementation:\n", + "1.039967668952797\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "1.0399676689527975\n", + "MSE values for own Ridge implementation\n", + "7.571105947979326e-05\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "7.57110594797945e-05\n", + "Beta values for own Ridge implementation\n", + "[-0.05100875 -0.04063602 -0.02723445 -0.01713366 -0.0100706 -0.00517114\n", + " -0.00174276 0.00068734 0.00243186 0.00369758 0.00462287 0.0053018\n", + " 0.00579953 0.006162 0.00642221 0.00660427 0.00672607 0.0068011\n", + " 0.00683964]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[-0.05100875 -0.04063602 -0.02723445 -0.01713366 -0.0100706 -0.00517114\n", + " -0.00174276 0.00068734 0.00243186 0.00369758 0.00462287 0.0053018\n", + " 0.00579953 0.006162 0.00642221 0.00660427 0.00672607 0.0068011\n", + " 0.00683964]\n", + "Intercept from own implementation:\n", + "0.999955585168597\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "0.999955585168597\n", + "MSE values for own Ridge implementation\n", + "0.0007698473260556339\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "0.000769847326055633\n", + "Beta values for own Ridge implementation\n", + "[-0.00834567 -0.00803064 -0.00673407 -0.00554552 -0.00458878 -0.0038335\n", + " -0.00323332 -0.00274989 -0.0023548 -0.00202756 -0.00175331 -0.00152117\n", + " -0.001323 -0.0011526 -0.00100519 -0.00087697 -0.00076495 -0.00066668\n", + " -0.00058016]\n", + "Beta values for Scikit-Learn Ridge implementation\n", + "[-0.00834567 -0.00803064 -0.00673407 -0.00554552 -0.00458878 -0.0038335\n", + " -0.00323332 -0.00274989 -0.0023548 -0.00202756 -0.00175331 -0.00152117\n", + " -0.001323 -0.0011526 -0.00100519 -0.00087697 -0.00076495 -0.00066668\n", + " -0.00058016]\n", + "Intercept from own implementation:\n", + "0.9637117593816477\n", + "Intercept from Scikit-Learn Ridge implementation\n", + "0.9637117593816477\n", + "MSE values for own Ridge implementation\n", + "0.0023813163025848865\n", + "MSE values for Scikit-Learn Ridge implementation\n", + "0.002381316302584886\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0tUlEQVR4nO3dd3xUVdrA8d+TDoQOIp2ICIYWpIkURUBRIthYiUoRfLHAsor6CjZEdC0oYAHWBjaWIliiUgTFF1dUitKLFFFBpISWEFLnef+YSzaEFAgzuZPk+X4+89m595x77nNGdp6ce+6cK6qKMcYYc66C3A7AGGNMyWAJxRhjjE9YQjHGGOMTllCMMcb4hCUUY4wxPmEJxRhjjE/4NaGISE8R2Soi20VkVC7l4SIy2yn/UUQaZCsb7ezfKiJXO/vqishSEdkkIhtF5B/Z6j8pIntEZI3zutaffTPGGHMq8dfvUEQkGPgF6AHsBlYCcaq6KVude4EWqnq3iPQDblDVW0QkGpgJtANqAUuAi4DzgJqq+pOIlAdWA9er6iYReRJIUtUX/dIhY4wx+fLnCKUdsF1Vd6pqGjAL6JOjTh/gXef9XKCbiIizf5aqpqrqr8B2oJ2q7lXVnwBUNRHYDNT2Yx+MMcacoRA/tl0b+CPb9m6gfV51VDVDRI4CVZ39P+Q49pTE4VweawX8mG33cBEZAKwCHlDVwzmDEpGhwFCAcuXKtW7SpMlZd8wYY0qz1atXH1TV6jn3+zOh+I2IRALzgPtU9ZizeyowDlDnf18CBuc8VlXfAN4AaNOmja5atapIYjbGmJJCRH7Lbb8/L3ntAepm267j7Mu1joiEABWBhPyOFZFQvMlkhqp+dLKCqu5T1UxV9QBv4r3kZowxpoj4M6GsBBqJSJSIhAH9gPgcdeKBgc77m4Gv1XuXQDzQz7kLLApoBKxw5lfeBjar6oTsDYlIzWybNwAbfN4jY4wxefLbJS9nTmQ4sAgIBqap6kYReQpYparxeJPD+yKyHTiEN+ng1JsDbAIygGGqmikinYD+wHoRWeOc6hFVnQ+8ICIxeC957QLu8lffjDHGnM5vtw0XB7nNoaSnp7N7925SUlJcisoEooiICOrUqUNoaKjboRjjOhFZraptcu4vlpPy/rR7927Kly9PgwYN8F5hM6WdqpKQkMDu3buJiopyOxxjApYtvZJDSkoKVatWtWRisogIVatWtVGrMQWwhJILSyYmJ/s3YUzBLKEYY0wpkpHp8VvbllACkIhw++23Z21nZGRQvXp1YmNjAdi3bx+xsbG0bNmS6Ohorr3Wuw7mrl27KFOmDDExMVmv9957r0hjzx5DdHQ0AwYMID09HYBVq1YxYsSIXI9r0KABBw8ePKdzT58+PavfYWFhNG/enJiYGEaNOm1d0jxNmjSJ5OTkc4rDmED2697D3DP1A3YfOFZw5bOlqqX21bp1a81p06ZNp+0rauXKldOWLVtqcnKyqqrOnz9fW7Zsqb169VJV1aFDh+qkSZOy6q9du1ZVVX/99Vdt2rRp0QecTfYYMjIytGvXrvrBBx8UeFz9+vX1wIEDPoujsO3ld1wg/NswJhDg/enHad+pNkIJUNdeey1ffPEFADNnziQuLi6rbO/evdSpUydru0WLFmfV9ldffUWrVq1o3rw5gwcPJjU1lZUrV3LjjTcC8Omnn1KmTBnS0tJISUnhggsuAOCKK67g4Ycfpl27dlx00UV8++23+Z4nODiYdu3asWePd4GEb775JmuUlZCQwFVXXUXTpk2588470Wy3r48bN47GjRvTqVMn4uLiePFF7wLSO3bsoGfPnrRu3ZrOnTuzZcuWM+rv+PHjadu2LS1atGDMmDEAHD9+nF69etGyZUuaNWvG7NmzeeWVV/jzzz/p2rUrXbt2PYtP1Jji4fX5y5n6xXd+a98SSgGuuOL015Qp3rLk5NzL33nHW37w4OllZ6pfv37MmjWLlJQU1q1bR/v2/11Xc9iwYQwZMoSuXbvyzDPP8Oeff2aV7dix45RLXjm/9FNSUhg0aBCzZ89m/fr1ZGRkMHXqVFq1asWaNWsA+Pbbb2nWrBkrV67kxx9/POXcGRkZrFixgkmTJjF27Nh8+5CSksKPP/5Iz549TysbO3YsnTp1YuPGjdxwww38/vvvAKxcuZJ58+axdu1aFixYQPbfCQ0dOpRXX32V1atX8+KLL3LvvfcW+Dl++eWXbNu2jRUrVrBmzRpWr17NsmXLWLhwIbVq1WLt2rVs2LCBnj17MmLECGrVqsXSpUtZunRpgW0bU9w8+OWDjFjan7T0TL+0b79DCVAtWrRg165dzJw5M2uO5KSrr76anTt3snDhQhYsWECrVq3YsMG70kzDhg2zEkNutm7dSlRUFBdddBEAAwcOZPLkydx33300bNiQzZs3s2LFCkaOHMmyZcvIzMykc+fOWcefHMW0bt2aXbt25XqOk0nt119/pVevXrmOoJYtW8ZHH3mXYuvVqxeVK1cG4LvvvqNPnz5EREQQERHBddddB0BSUhLLly+nb9++WW2kpqbm9xEC3oTy5Zdf0qpVq6x2tm3bRufOnXnggQd4+OGHiY2NPaWPxpREU7/4jqTK33Nz2VcJCw32yzksoRTgm2/yLitbNv/yatXyLy9I7969efDBB/nmm29ISEg4paxKlSrceuut3HrrrcTGxrJs2TJat25d+JMBXbp0YcGCBYSGhtK9e3cGDRpEZmYm48ePz6oTHh4OeC9nZWRk5NrOyaR28OBBOnbsSHx8PL179z6n2DweD5UqVco3WeZGVRk9ejR33XX6Sjw//fQT8+fP57HHHqNbt2488cQT5xSjMYFs3JIXkfAqTB5xh9/OYZe8AtjgwYMZM2YMzZs3P2X/119/nXUnUmJiIjt27KBevXpn1Gbjxo3ZtWsX27dvB+D999/n8ssvB6Bz585MmjSJDh06UL16dRISEti6dSvNmjUrVPzVqlXjueee49lnnz2trEuXLvz73/8GYMGCBRw+7H10TceOHfnss89ISUkhKSmJzz//HIAKFSoQFRXFhx9+CHgTxdq1awuM4eqrr2batGkkJSUBsGfPHvbv38+ff/5J2bJluf3223nooYf46aefAChfvjyJiYmF6q8xgWrRql/YW/FTOobfy3mVy/ntPDZCCWB16tTJ9Tbb1atXM3z4cEJCQvB4PNx55520bduWXbt2ZV1uOmnw4MGntBEREcH06dPp27cvGRkZtG3blrvvvhuA9u3bs2/fPrp06QJ4L7v99ddf5/Sjvuuvv54nn3zytLmcMWPGEBcXR9OmTbnsssuyEmLbtm3p3bs3LVq0oEaNGjRv3pyKFSsCMGPGDO655x6efvpp0tPT6devHy1btsz3/FdddRWbN2+mQ4cOAERGRvLBBx+wfft2HnroIYKCgggNDWXq1KmAd56mZ8+eWXMpxpQE3235haDjtZh673C/nscWh8yxOOTmzZu5+OKLXYrIgHeeIzIykuTkZLp06cIbb7zBJZdc4nZY9m/DFGtp6Zk+mzvJa3FIu+RlAs7QoUOJiYnhkksu4aabbgqIZGJMcTV/xRafJpP82CUvE3BOzq0YY87NwaPJxH7UiSZzb2bTC//y+/ksoRhjTAn197feRcskcO8ltxXJ+SyhGGNMCZSWnsm8PydQztOOe3t1KpJzWkIxxpgS6PEZ8aRX2M7f63xIUFDRPH7BEooxxpRAM9fNJiQoimcH3FBk57S7vAKQv5ev//zzz2nVqlXW8a+//joATz75ZNZCjHlp0KABzZs3p0WLFlx++eX89ttvWWWXXXZZrscMGjSIuXPnnt2HkMP69euz+lSlShWioqKIiYmhe/fuZ9zGJ598wqZNm84pDmOKi+3Pz2Bx/yVFcnfXSTZCCUDlypVjw4YNnDhxgjJlyrB48WJq166dVf7EE0/Qo0cP/vGPfwCwbt26rLKC1vJKT09n6NChrFixgjp16pCamprnmlx5Wbp0KdWqVWPMmDE8/fTTvPnmmwAsX778rNo5G82bN8/q16BBg4iNjeXmm28+qzY++eQTYmNjiY6O9kOExgSO5JR0ykaEckXLC4r0vDZCCVD+Wr4+MTGRjIwMqlatCnjX5mrcuHGhYuzQoUPW0vTg/RU6eJdFGT58OI0bN6Z79+7s378/q878+fNp0qQJrVu3ZsSIEVmjruPHjzN48GDatWtHq1at+PTTT88ohi+//JIOHTpwySWX0Ldv36wlVkaNGkV0dDQtWrTgwQcfZPny5cTHx/PQQw8RExPDjh07CtVnYwLdolW/UH5MXSZ8XPQrPdgIpQBXvHPFafv+1vRv3Nv2XpLTk7l2xrWnlQ+KGcSgmEEcTD7IzXNO/Sv6m0HfnNF5+/Xrx1NPPUVsbCzr1q1j8ODBWcuXDBs2jFtuuYXXXnuN7t27c8cdd1CrVi2A05ZeefXVV09ZSbdKlSr07t2b+vXr061bN2JjY4mLiyMo6Oz/tli4cCHXX3/9afs//vhjtm7dyqZNm9i3bx/R0dEMHjyYlJQU7rrrLpYtW0ZUVNQpSfKZZ57hyiuvZNq0aRw5coR27drRvXt3ypXLe92hgwcP8vTTT7NkyRLKlSvH888/z4QJExg2bBgff/wxW7ZsQUQ4cuQIlSpVonfv3oUa2RhTnIz8cCKesCNcFVP0I3FLKAHKX8vXA7z11lusX7+eJUuW8OKLL7J48WLeOfkQlzPQtWtXDh06RGRkJOPGjTutfNmyZcTFxREcHEytWrW48sorAdiyZQsXXHABUVFRAMTFxfHGG28A3pFGfHx81hxOSkoKv//+e75Lnfzwww9s2rSJjh07ApCWlkaHDh2oWLEiERERDBkyhNjY2KxRkDEl3cZd+9kU+g5NUgfQLKpGkZ/fEkoB8htRlA0tm295tbLVznhEkht/Ll/fvHlzmjdvTv/+/YmKijqrhLJ06VIqVarEbbfdxpgxY5gwYcIZH5sXVWXevHlndflNVenRowczZ848rWzFihV89dVXzJ07l9dee42vv/76nGM0JtANe2cKhKbwUu+Rrpzf5lACmD+Wr09KSuKbbA9pWbNmDfXr1z/r2EJCQpg0aRLvvfcehw4dOqWsS5cuzJ49m8zMTPbu3Zu1am/jxo3ZuXNn1k0As2fPzjrm6quv5tVXX816FPDPP/9cYAyXXnop3333XdZS/MePH+eXX34hKSmJo0ePcu211zJx4sSsZe5taXpTkh1JSmFZymRqHLmOa9s1cSUGSygBLL/l69u0aUOLFi3o0KFD1vL1cPojgF955ZVTjlVVXnjhBRo3bkxMTAxjxow5ZXTy9NNPU6dOnaxXfmrWrElcXByTJ08+Zf8NN9xAo0aNiI6OZsCAAVlLx5cpU4YpU6ZkPRe+fPnyWUvTP/7446Snp9OiRQuaNm3K448/XuDnU716dd555x3i4uKyPostW7aQmJhIbGwsLVq0oFOnTlkjqH79+jF+/HhatWplk/KmxKlQNpyXL5vHKzecfhm6qNjy9bZ8fZE6uTS9qjJs2DAaNWrE/fff73ZYZ8T+bRjjZcvXm4Dw5ptvEhMTQ9OmTTl69Giuj+Y1xpydx9//jBajhrP7wDFX47BJeVOk7r///mIzIjGmuJi0+llSQv7ivEovuxqHJZRcqOo5PfbWlDyl+dKwCWyvz19OUuXvubnsq0W6zEpu7JJXDhERESQkJNgXiMmiqiQkJBAREeF2KMacZuzi8ciJKky+8w63Q7ERSk516tRh9+7dHDhwwO1QTACJiIgo8K43Y4raolW/sLfip3TSRzmvct6rShQVvyYUEekJvAwEA2+p6nM5ysOB94DWQAJwi6rucspGA0OATGCEqi4SkbpO/RqAAm+o6stO/SrAbKABsAv4m6oePtuYQ0NDs37JbYwxgSyyTDiNjw9h6r3D3Q4F8OMlLxEJBiYD1wDRQJyI5FxcZghwWFUvBCYCzzvHRgP9gKZAT2CK014G8ICqRgOXAsOytTkK+EpVGwFfOdvGGFNidWxany3j33RlmZXc+HMOpR2wXVV3qmoaMAvok6NOH+Bd5/1coJt4Z8P7ALNUNVVVfwW2A+1Uda+q/gSgqonAZqB2Lm29C1zvn24ZY4z7Hnh7LtO/XOF2GKfwZ0KpDfyRbXs3//3yP62OqmYAR4GqZ3KsiDQAWgE/OrtqqOpe5/1feC+LnUZEhorIKhFZZfMkxpji6ODRZCZuu5vRC552O5RTFMu7vEQkEpgH3Keqp/2SR723aOV6m5aqvqGqbVS1TfXq1f0cqTHG+N7f33oXLZPAo10fdDuUU/gzoewB6mbbruPsy7WOiIQAFfFOzud5rIiE4k0mM1T1o2x19olITadOTWA/xhhTwqSlZzLvzwmUO9KOYbGdCz6gCPkzoawEGolIlIiE4Z1kj89RJx4Y6Ly/GfjaGV3EA/1EJFxEooBGwApnfuVtYLOq5lwzPXtbA4Eze+SfMcYUI4998CnpFbYztNmDBAUF1g+w/ZZQnDmR4cAivJPnc1R1o4g8JSK9nWpvA1VFZDswEufOLFXdCMwBNgELgWGqmgl0BPoDV4rIGud18ulTzwE9RGQb0N3ZNsaYEmXfsUOUPdKa5wbc6HYop7HVhnOsNmyMMYHO41FXRye22rAxxhRz7y1ZRUamJ+AudZ1kCcUYY4qBRat+YeB/2nHDCxPdDiVPllCMMaYYGPnhRMgM5dl+t7sdSp4soRhjTIDbuGs/m0LfoXHKgIBZZiU3ttqwMcYEuGHvTIHQFF68bqTboeTLEooxxgQwj0dZefRzanAdse0vdjucfFlCMcaYABYUJBx49nt+/eusn8ZR5CyhGGNMgEpLzyQlLYMK5cJp2uA8t8MpkE3KG2NMgHp8RjyVx0axYOVWt0M5I5ZQjDEmQE1ZM54gTwRdWzZ0O5QzYpe8jDEmAL0+fzlJlb/nprKvEBFWPL6qbYRijDEBaOziF5GUyky5c7DboZwxSyjGGBNglq37lb0VP+Gy0Hs5r3I5t8M5Y8VjHGWMMaVIp2YNeGHbYq65pJnboZwVSyjGGBNggoKEh27q5nYYZ80ueRljTADp9c/xtH5kJBmZHrdDOWuWUIwxJkAcPJrMgmMvsPvENkKCi9/Xs13yMsaYADHi7ffQMgd57LKH3A6lUCyhGGNMAEhLz2Tunpco52nLsNjObodTKJZQjDEmADw+I570CtsZVnt2wD7ityDF7yKdMcaUQJ0vbkKrlPt4dsCNbodSaDZCMcaYABDb/mJi2wfu8+LPhI1QjDHGZTe+8DJzlq11O4xzZgnFGGNctHj1Nj5Ovp9XlnzodijnzBKKMca46L45EyAzlCmDhrsdyjmzhGKMMS7Z/PsBNoW+w0Up/Wlxwfluh3PObFLeGGNccu/0KRCawkvXPeB2KD5hIxRjjHFJSFAI9Y/dSmz7i90OxSdshGKMMS5Z/PijbofgUzZCMcaYIpaWnsmEj5fi8ajbofiUJRRjjClij8+I54F1V/Lkv79wOxSfsktexhhTxKaueZGQkAY88reebofiU34doYhITxHZKiLbRWRULuXhIjLbKf9RRBpkKxvt7N8qIldn2z9NRPaLyIYcbT0pIntEZI3zutaffTPGmMJ4ff5yEisvp0+NkUSElay/6f2WUEQkGJgMXANEA3EiEp2j2hDgsKpeCEwEnneOjQb6AU2BnsAUpz2Ad5x9uZmoqjHOa74v+2OMMb4wdvGLSEplXhtyh9uh+Jw/RyjtgO2qulNV04BZQJ8cdfoA7zrv5wLdRESc/bNUNVVVfwW2O+2hqsuAQ36M2xhj/OKvQ0nsD1lJh9B7OL9KpNvh+Jw/x1u1gT+ybe8G2udVR1UzROQoUNXZ/0OOY2ufwTmHi8gAYBXwgKoezllBRIYCQwHq1at3Zj0xxhgfOL9KJEee2kHSiTS3Q/GLknSX11SgIRAD7AVeyq2Sqr6hqm1UtU316tWLMDxjTGl28GgyySnpRJYJK5GjE/BvQtkD1M22XcfZl2sdEQkBKgIJZ3jsKVR1n6pmqqoHeBPnEpkxxgSCvpPGU+HxC/kzIdHtUPzGnwllJdBIRKJEJAzvJHt8jjrxwEDn/c3A16qqzv5+zl1gUUAjYEV+JxORmtk2bwA25FXXGGOK0qFjJ/i/E69RNaMFtaqWdzscv/HbHIozJzIcWAQEA9NUdaOIPAWsUtV44G3gfRHZjneivZ9z7EYRmQNsAjKAYaqaCSAiM4ErgGoishsYo6pvAy+ISAygwC7gLn/1zRhjzsbwt95FyxzkkQ4Puh2KX4l3QFA6tWnTRletWuV2GMaYEiwtPZPIURcT5qnEsZd+JChI3A7pnInIalVtk3N/yfpVjTHGBJgX5i0mvcI2htWeXSKSSX4soRhjjB898rerCZn3Jff16ep2KH5nCcUYY/woKEgY1beH22EUiZL0OxRjjAkoFz00hKvG/dPtMIqMJRRjjPGDxau3sa3cdI6nH3c7lCJjCcUYY/zg/jkTITOUyQOHux1KkbGEYowxPrb59wNsDJ3ORSn9iWlYs+ADSgiblDfGGB+7d/oUCE3hpesecDuUImUJxRhjfGzo5b0pu7wsse0vdjuUImUJxRhjfCzuilbEXdHK7TCKnM2hGGOMj6SlZ9L+sYdZsHKr26G4whKKMcb4yJh/f8aK0Bf4cu1at0NxRb4JRURuz/a+Y46y0nMvnDHGnIHJP79ISGIDnh1wo9uhuKKgEcrIbO9fzVE22MexGGNMsfXGgu9JrPwdfWqMJCKsdE5PF5RQJI/3uW0bY0ypNXbxi0hKZV4bcofbobimoDSqebzPbdsYY0olj0epEV6P6LIPlNjnxZ+JghJKExFZh3c00tB5j7N9gV8jM8aYYiIoSPjp2Yluh+G6ghJK6fpVjjHGnKWtfxxkweoNjOh9eYl/gFZB8p1DUdXfsr+AJOASoJqzbYwxpdo906Zw/9quLF27w+1QXFfQbcOfi0gz531NYAPeu7veF5H7/B+eMcYErkPHTvBN8mucdySWbq0udDsc1xV0l1eUqm5w3t8BLFbV64D22G3DxphS7u9vv4eWPcAjVzzodigBoaCEkp7tfTdgPoCqJgIefwVljDGBLiPTw4e7X6LskTb8/boubocTEAqalP9DRP4O7MY7d7IQQETKAKF+js0YYwLW0jU7yAg5yrBGT5f6yfiTCkooQ4CngO7ALap6xNl/KTDdj3EZY0xA69G6EYca/1ZqfxWfm3w/CVXdD9ydy/6lwFJ/BWWMMYFs6x8HqV+jEpUiI9wOJaDkm1BEJD6/clXt7dtwjDEm8F058S4O6y6SXlpll7uyKWis1gH4A5gJ/Iit32WMKeW++nk7f1b4mMs8oy2Z5FBQQjkf6AHEAbcCXwAzVXWjvwMzxphAdN+siRAWytQ7/u52KAGnoF/KZ6rqQlUdiHcifjvwjT0LxRhTGm394yAbQqdzUUp/WlxwvtvhBJwCb08QkXCgF95RSgPgFeBj/4ZljDGBZ9S/P4DQE4yPHVlw5VKooEn594BmeH/QODbbr+aNMabU+fDBEby9qB29L412O5SAJKp5P9ZERDzAcWcze0UBVFUr+DE2v2vTpo2uWrXK7TCMMcWAx6M2Ce8QkdWq2ibn/oLmUIJUtbzzqpDtVb64JxNjjDlTGZkeKo7syMCX33I7lIBW0Fpe50REeorIVhHZLiKjcikPF5HZTvmPItIgW9loZ/9WEbk62/5pIrJfRDbkaKuKiCwWkW3O/1b2Z9+MMaXHEzM+I6ny91QpZ39H58dvCUVEgoHJwDVANBAnIjkvPA4BDqvqhcBE4Hnn2GigH9AU6AlMcdoDeMfZl9Mo4CtVbQR85WwbY8w5e+2n8YQkNuDZATe6HUpA8+cIpR2wXVV3qmoaMAvok6NOH+Bd5/1coJuIiLN/lqqmquqveG9XbgegqsuAQ7mcL3tb7wLX+7AvxphS6o0F35NY+Tt6n3e/rdtVAH8mlNp4f2V/0m5nX651VDUDOApUPcNjc6qhqnud938BNXKrJCJDRWSViKw6cODAmfTDGFOKPbX4JSSlEpPvtEdAFcSvcyhuUe+ta7nevqaqb6hqG1VtU7169SKOzBhT3DzR4wHuqT+Z86tEuh1KwPPn+G0PUDfbdh1nX251dotICFARSDjDY3PaJyI1VXWv87ji/ecSvDHGAAy9pgPeZQ1NQfw5QlkJNBKRKBEJwzvJnnP14nhgoPP+ZuBrZ3QRD/Rz7gKLAhoBKwo4X/a2BgKf+qAPxphSausfB2n28L18v+l3t0MpNvyWUJw5keHAImAzMEdVN4rIUyJyctn7t4GqIrIdGIlzZ5az+OQcYBPep0QOU9VMABGZCXwPNBaR3SIyxGnrOaCHiGzD+0Cw5/zVN2NMyXfPtClsLDuVhMTjBVc2QAG/lC/p7JfyxpjcHDp2gmrP1Kd6Wnv2TfzM7XACTl6/lLd74IwxJocR095Hyx7gkcsedDuUYsUSijHGZJOR6WHOHy9R1tOGv1/Xxe1wihVLKMYYk83+I8dpFHIlsRdfZYtBniVLKMYYk02tquXZ+PxUt8MolkrkDxuNMaYwPl2+kalffOd2GMWWJRRjjHEMmzuGYf+5joNHk90OpViyhGKMMcBXP29nT4WP6BB6D9UqlnU7nGLJEooxxgD3zZoImaFMHjjc7VCKLUsoxphSb+sfB9kQOp2LUvoT07Cm2+EUW5ZQjDGl3oLVG5CMcoy/aaTboRRrtvSKLb1ijAGSTqQRWSbM7TCKhbyWXrERijGmVPt+0++kpWdaMvEBSyjGmFLr2PFULn+7Bw0fjnM7lBLBEooxplTaf/g4UY9eR3qFXxjQ6na3wykRbOkVY0yp8/v+ozT/ZyzHKi7njirTeKZ/74IPMgWyhGKMKXUufXYQxyr8wP31ZjFhSF+3wykxLKEYY0qd6bc9z4/bhvJE3DVuh1Ki2ByKMaZU+M+GXXR7ahwej3J1m4ssmfiBJRRjTIk3f8UWLn+3E0tPTOD7zb+7HU6JZQnFGFOizVm2lth5XdCgdOZc+390bFrf7ZBKLEsoxpgS662FP9BvwRUEecL54uZl3Ny5hdshlWg2KW+MKbF2JyQQmlqTr4cssJFJEbC1vGwtL2NKnDU79matGpySlkFEmP3t7Eu2lpcxplS4/605tJp+AU/PWghgyaQI2SdtjCkx7nhlGu8k/A8VkjoyqNtlbodT6lhCMcaUCDe98AofnfgHVY9exZaxH9tjfF1gl7yMMcXelM//w0cn/kHNIzew85l4SyYusRGKMabYu/vajmzdO4vnR99kcyYushGKMaZYysj00P6xh5n3n/UEBQkv/88tlkxcZgnFGFPspKRl0PjhQawIfYF/LY13OxzjsHRujClWjh1P5eLHb+XPih/RTcax6NFH3A7JOCyhGGOKjYNHk2ny5I0kVFrE9RET+fjh+9wOyWTj10teItJTRLaKyHYRGZVLebiIzHbKfxSRBtnKRjv7t4rI1QW1KSLviMivIrLGecX4s2/GmKIXFhKMIAys/JYlkwDktxGKiAQDk4EewG5gpYjEq+qmbNWGAIdV9UIR6Qc8D9wiItFAP6ApUAtYIiIXOcfk1+ZDqjrXX30yxrhj2+4EgoKEhrWqsO+l+QQFidshmVz4c4TSDtiuqjtVNQ2YBfTJUacP8K7zfi7QTUTE2T9LVVNV9Vdgu9PembRpjClB1uzYS/MJl9P6hRvweNSSSQDzZ0KpDfyRbXu3sy/XOqqaARwFquZzbEFtPiMi60RkooiE+6ITxhj3fLfxN9pN6UJqmV083mWMJZMAV5JuGx4NNAHaAlWAh3OrJCJDRWSViKw6cOBAUcZnjDkLi1b9wuXvdCIj7CBvdl7CAzde6XZIpgD+TCh7gLrZtus4+3KtIyIhQEUgIZ9j82xTVfeqVyowHe/lsdOo6huq2kZV21SvXr2QXTPG+JPHo9z8wSA8QanM7LmUO3te6nZI5gz487bhlUAjEYnC+6XfD7g1R514YCDwPXAz8LWqqojEA/8WkQl4J+UbASsAyatNEampqnudOZjrgQ1+7Jsxxo+CgoQvBs8gKSWVa9s1cTscc4b8llBUNUNEhgOLgGBgmqpuFJGngFWqGg+8DbwvItuBQ3gTBE69OcAmIAMYpqqZALm16ZxyhohUx5t01gB3+6tvxhj/mPTJN0z7YR4/PfMyXVpEuR2OOUv2xEZ7YqMxAWHsv+fz5KabCEuOYtOD/6FhrSpuh2TyYE9sNMYErJFvf8iTW66nzPFo1t2/zJJJMWVLrxhjXDV08nu8uf8Oyh/rwIZHvqDeeRXdDskUko1QjDGuurBGLc47di3bxyyyZFLMWUIxxrhixtc/AfC/N3dn38TPOK9yOZcjMufKEooxpkh5PEqHx0dz+7eteTV+mdvhGB+yORRjTJHJyPRwyaMjWF9mMhcfv4t7enVyOyTjQzZCMcYUiZS0DJo8PJj1ZSbTOu0BNjw3lZBg+woqSey/pjGmSIybPZ8d5d/lSnmKFePG20KPJZBd8jLGFIln+vem/oLlDL2mg9uhGD+xEYoxxm92HzjG+ff35oOvVgNYMinhLKEYY/xix5+HaPLP7uwrv4Cfdu10OxxTBOySlzHG59bt/It2r/UgNXIbjzX6iHH9r3M7JFMELKEYY3xq3c6/aDO5C+ll9vBCzBc8dFM3t0MyRcQueRljfOrCWlWpQ3v+1XGxJZNSxkYoxhif+HT5Ri6sVZ2mDc5j50vvux2OcYGNUIwx5+zdxSu5Ib4LXV++w+1QjIssoRhjzsmr8csYtLQbwRkVmTPoVbfDMS6yhGKMKbSnZy1kxIqehKXU4fuh33JFywvcDsm4yOZQjDGFkpaeyTMr/pcy2oTV9y/i4nrV3Q7JuMwSijHmrHk8SlhoMN/evYDqFctRv0Ylt0MyAcAueRljzkq/l6ZwwUO3k5aeSZuLalsyMVksoRhjztg1zzzP7KRhpHqSSEnLcDscE2AsoRhjCuTxKB2feJSFGaOofyyOHc/OpUK5cLfDMgHGEooxpkCdxjzK8uB/0uT4//DLc+9TNiLU7ZBMALJJeWNMgW5tdw26QvnuuX/ag7FMnmyEYozJVdKJNJ744HMAhl/Xme/HPWvJxOTLEoox5jSHjp2g4SM3Mm7Hdcz7z3q3wzHFhF3yMsac4s+ERJqO68ORSt9wW8XXualTc7dDMsWEJRRjTJZf9x6m+QvXcLziKu45/32m3H2b2yGZYsQSijEBKC09k2PJqSQmp3I0OYUyYaE0rlsNgA++Ws2R5GSSU1NJTk3leGoKzerWY2CPtng8yk3jXyYlI5XUjFRSM72vHhd14tmB13PwaDKtxg4iQ1NJJ5UMTSFTUulT/w4+uG8oE+IXcjxyDaMbzuOfA/q4/CmY4sYSijE5JJ1IY8/BYxw4msSR4ydIPJFCkAi3XB4DwOvzl7N17x5OpKVyIj2VlPRUKpetwNR7bgeg7/jX2HFoB2meVNI8KaR5UqkX2ZBlY58C4IIH+rPf8wuZpOKRVDxBqdSjEztefA+AkAcbkFn+t1NiqnO0L39MmANA/yXdIOLoKeWNNt/BwB5tCQoSPkl+EIIynZIQIJyQbSHA9YQEB7Gf9QQRTrDzCtPylA2NAODVu+Lot/EyOjat7/sP1pR4llBMsZaSlsG+w0nsP5JE60a1CQoS5q/Ywndbf+HI8SSOpiRx9EQiJ9JPsOSJxwCIe2kq3+xeSKomkSZJZEgSoqGcmLgGgIse6cfeSh+fcp7gpHrccrn3S/7RL8eSUPnLU8rD90QzFW9CWbJnHkfKrkYknCAiEMLxJGZm1Q0NCifCU5lQwgklglDCaVYtJqv8mqp/53jaccJDwokICadMaAStmjXKKn/2kg8JEqFseDjlIyIoFxHOBedXyyrfefcBypcNp0LZcMJCg0+Js1JkBKkTNuf7mVoyMYVlCcUUub8OJbH+170cTEzi4LFEDiUlceh4Ig9dfw11qlfg9fnLmb48nuSMJI5nJJLiSSLFk8gPD86kUZ2q9HluAp8deRYNSYLQlKx29ww/Rq2q5Xni07dYHfbSaedNSx9NWGgwe47t5Qi/EUokEVQmgnqUD62SVe+OSwaxZe8VVIiIpHxEWcqFR1CtfIWs8tkDJ3MsOYXIiHAiy4RTvkw4lSLLZJUfnrQ03/5vHf9WvuWfjX4g3/JRfXvkWx5Vs3K+5cb4i18Tioj0BF4GgoG3VPW5HOXhwHtAayABuEVVdzllo4EhQCYwQlUX5demiEQBs4CqwGqgv6qm+bN/JUlauvcv6LDQYJJT0tn0+35S0tI5npJGcqr3dWmTBjSsVYU1O/YyeeGXHElOJDE1iaS0JI6nJzG2z1B6XxrNq/HLeOzrx0iXRDKCksgMScQTksSbnZdwZ89LeWTGh0w/NPi0GDpvWcvN1Vvw1aaf+VEmIpQniEhCiCREIkk8kQpAq7qN2XG4L2U1kvJSnvJhkVSIiKRsuPfX25PihrNrXxxVK0RSo1J5zqsUyXmVymX9te699PRUnp/FM/175/tZdWt1YWE+YmNKPL8lFBEJBiYDPYDdwEoRiVfVTdmqDQEOq+qFItIPeB64RUSigX5AU6AWsERELnKOyavN54GJqjpLRP7ltD3VX/0rSNKJNJJOpHEi1fulfDwljcgy4TRtcB4AM77+iaSUFE6kpXEiLY3U9HQurlObWy6PweNRBr82jbSMdFIz0kjNTCMtI40rm7RlVN8eHElKoeszo0j3pJHhSSdd08jQNG68+AYm3vk3tu1OoN2LfckkDY+k4SEdj6QR1+ABpo8YzFc/b6fH7EvRoDQIdl6i3FbxdT64byhz/7OWgd+1Pa1Pw7f9m1fviuOb9Vt56+Cg/xZIEFCedbt60vvSaEKCgwkihArUJYJIymh5ymkkF5zvfV7G7Z27EPHDB1QqG0nlcpFUiYykesXydGnmfTjTrJH3MidoWJ6f7ZO39eLJ23rlWd6pWQM6NWtw1v/NjDHnxp8jlHbAdlXdCSAis4A+QPaE0gd40nk/F3hNRMTZP0tVU4FfRWS70x65tSkim4ErgVudOu867fotoZx3/7UcDt6KShqeoDQISqNm6pXsmTAPgIpPXIAncs8px9Q9+jd+nzAbgNsXX3n6xOqawdxy+dsEBQnvJgyFIM8p5Qd++gej+vbA41HWyHREwhAJI0jDEEL57dClAIQEB+EhgxDCCaY8IYQRTBi1Knuvs9eqWpFm9COUMEI1jFANJSwojJ7NWwPQvnEDbtvwOhEhYYSHhBEeGkpEaBh92l7ijb1rO2KidlCjcnlqVIqkUmTEKb+gvqdXR+7p9XWen92VMQ25MqZhnuX2a2xjiid/JpTawB/ZtncD7fOqo6oZInIU7yWr2sAPOY6t7bzPrc2qwBFVzcil/ilEZCgwFKBevXpn16Ns6pdtSmRKVUIkjFDCCNFQmtVtmlV+U43RHE9L9n4hB4cRHhpGTMv/fok+HTObIBHCQ0MpExZG2bAwLqx1Xlb58lt+pVxEGOUiwigTHkpkRFjWgnxVKpRB/3lqMsouqmZljk5almf5xfWqs+651/Isb1y3Gh/cNzTP8moVy9qjXo0xpyl1k/Kq+gbwBkCbNm20sO2sfGZ8vuVzHsz7kg3Ao7dcnW95h+jCJztjjHGDP9fy2gPUzbZdx9mXax0RCQEq4p2cz+vYvPYnAJWcNvI6lzHGGD/yZ0JZCTQSkSgRCcM7yR6fo048MNB5fzPwtaqqs7+fiIQ7d281Albk1aZzzFKnDZw2P/Vj34wxxuTgt0tezpzIcGAR3lt8p6nqRhF5ClilqvHA28D7zqT7IbwJAqfeHLwT+BnAMFXNBMitTeeUDwOzRORp4GenbWOMMUVEvH/cl05t2rTRVatWuR2GMcYUKyKyWlXb5Nxvz0MxxhjjE5ZQjDHG+IQlFGOMMT5hCcUYY4xPlOpJeRE5APxWYMXcVQMO+jAcN1lfAk9J6QdYXwLVufSlvqpWz7mzVCeUcyEiq3K7y6E4sr4EnpLSD7C+BCp/9MUueRljjPEJSyjGGGN8whJK4b3hdgA+ZH0JPCWlH2B9CVQ+74vNoRhjjPEJG6EYY4zxCUsoxhhjfMISig+IyAMioiJSze1YCktExonIOhFZIyJfikgtt2MqDBEZLyJbnL58LCKV3I6psESkr4hsFBGPiBTLW1VFpKeIbBWR7SIyyu14CktEponIfhHZ4HYs50JE6orIUhHZ5Pzb+ocv27eEco5EpC5wFfC727Gco/Gq2kJVY4DPgSdcjqewFgPNVLUF8Asw2uV4zsUG4EYg7+c5BzARCQYmA9cA0UCciES7G1WhvQP0dDsIH8gAHlDVaOBSYJgv/5tYQjl3E4H/BYr13Q2qeizbZjmKaX9U9UtVzXA2f8D79M5iSVU3q+pWt+M4B+2A7aq6U1XTgFlAH5djKhRVXYb3mU3FmqruVdWfnPeJwGagtq/aL3XPlPclEekD7FHVtSLidjjnTESeAQYAR4GuLofjC4OB2W4HUYrVBv7Itr0baO9SLCYHEWkAtAJ+9FWbllAKICJLgPNzKXoUeATv5a5iIb++qOqnqvoo8KiIjAaGA2OKNMAzVFA/nDqP4h3ezyjK2M7WmfTFGF8TkUhgHnBfjqsT58QSSgFUtXtu+0WkORAFnByd1AF+EpF2qvpXEYZ4xvLqSy5mAPMJ0IRSUD9EZBAQC3TTAP+h1Vn8NymO9gB1s23XcfYZF4lIKN5kMkNVP/Jl25ZQCklV1wPnndwWkV1AG1UtliuRikgjVd3mbPYBtrgZT2GJSE+8c1qXq2qy2/GUciuBRiIShTeR9ANudTek0k28f/2+DWxW1Qm+bt8m5c1Jz4nIBhFZh/cynk9vJyxCrwHlgcXOLdD/cjugwhKRG0RkN9AB+EJEFrkd09lwbo4YDizCO/k7R1U3uhtV4YjITOB7oLGI7BaRIW7HVEgdgf7Alc7/P9aIyLW+atyWXjHGGOMTNkIxxhjjE5ZQjDHG+IQlFGOMMT5hCcUYY4xPWEIxxhjjE5ZQjHGISNI5HDvcWVH3lFWnxesVp2ydiFySraymiHzuvL/i5PtzJSLfnMnqxCKyq6AVskVkiYhU9kVcpuSzhGKMb3wHdAd+y7H/GqCR8xoKTM1WNhJ4s0iiK7z3gXvdDsIUD5ZQjMnBGVWMd37ouV5EbnH2B4nIFOd5K4tFZL6I3Aygqj+r6q5cmusDvKdePwCVRKSmU3YTsDCX87cTke9F5GcRWS4ijZ39g0TkE+fcu5xR0Uin3g8iUiVbM/2dH61tEJF2zvFVnWfdbBSRtwDJds5PRGS1UzY0WzvxQFxhP0tTulhCMeZ0NwIxQEu8o47xThK4EWiA99ke/fH+gr0gua24W9tZjuSwqqbmcswWoLOqtsL7XJp/Zitr5sTRFngGSHbqfY93peiTyjrPtrkXmObsGwP8R1WbAh8D9bLVH6yqrYE2wAgRqQqgqoeB8JPbxuTH1vIy5nSdgJmqmgnsE5H/w/sF3gn4UFU9wF8isvQczlETOJBHWUXgXRFphPe5NKHZypY6z7FIFJGjwGfO/vVAi2z1ZoL3OR4iUkG8T67sgjcZoapfiMjhbPVHiMgNzvu6eC/RJTjb+4Fa2baNyZWNUIzxr7xW3D0BRORxzDi8iaMZcF2OetlHNJ5s2x5O/QMx55pKea6xJCJX4B2JdVDVlsDPOc4Z4cRrTL4soRhzum+BW0QkWESq4/3LfgXeifebnLmUGsAVZ9BWPDDAmZe5FDiqqnvxPp64QR7HVOS/y7wPKmQfTs77dHLOeRTvo4RvdfZfA5y8e6si3stvySLSBO+jYXHqCd7ntewqZBymFLGEYszpPgbWAWuBr4H/dZ5xMw/vHMgm4APgJ7xPt0RERjgrA9cB1jmT3uB9rsxOYDveO7ruBVDV48AOEbkwl/O/ADwrIj9T+MvSKc7x/wJOrow7FugiIhvxXvr63dm/EAgRkc3Ac3gfnXxSa+CHbI9VNiZPttqwMWdBRCJVNcmZpF4BdCzsA9WcOYvWqvqYT4P0IRF5GYhX1a/cjsUEPpuUN+bsfO5McIcB487l6Zyq+nExuHtqgyUTc6ZshGKMMcYnbA7FGGOMT1hCMcYY4xOWUIwxxviEJRRjjDE+YQnFGGOMT/w/KpKrMndbP4YAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "import numpy as np\n", "import pandas as pd\n", @@ -1120,10 +1231,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -1236,10 +1344,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "X = np.zeros((n, L ** 2))\n", @@ -1303,10 +1408,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "X_train_own = np.concatenate(\n", @@ -1322,10 +1424,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n", @@ -1410,10 +1509,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n", @@ -1424,10 +1520,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "beta = ols_svd(X_train_own,y_train)" @@ -1443,10 +1536,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "J = beta[1:].reshape(L, L)" @@ -1462,10 +1552,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -1531,10 +1618,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -1640,10 +1724,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "X = np.zeros((n, L ** 2))\n", @@ -1673,10 +1754,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "clf = skl.LinearRegression().fit(X_train, y_train)" @@ -1692,10 +1770,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "J_sk = clf.coef_.reshape(L, L)" @@ -1711,10 +1786,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -1769,10 +1841,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "_lambda = 0.1\n", @@ -1823,10 +1892,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n", @@ -1860,10 +1926,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "lambdas = np.logspace(-4, 5, 10)\n", @@ -1926,10 +1989,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(20, 14))\n", @@ -2120,10 +2180,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# Common imports\n", @@ -2194,10 +2251,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "agegroupmean = np.array([0.1, 0.133, 0.250, 0.333, 0.462, 0.625, 0.765, 0.800])\n", @@ -2281,10 +2335,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "\"\"\"The sigmoid function (or the logistic curve) is a\n", @@ -2703,10 +2754,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -2749,10 +2797,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -2815,10 +2860,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)" @@ -2834,10 +2876,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "correlation_matrix = cancerpd.corr().round(1)" @@ -2860,10 +2899,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -3348,7 +3384,25 @@ ] } ], - "metadata": {}, + "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.8.5" + } + }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/doc/pub/week39/html/week39-bs.html b/doc/pub/week39/html/week39-bs.html index 4af09f5c3..746b595ba 100644 --- a/doc/pub/week39/html/week39-bs.html +++ b/doc/pub/week39/html/week39-bs.html @@ -48,6 +48,7 @@ Automatically generated HTML file from DocOnce source None, 'searching-for-optimal-regularization-parameters-lambda'), ('Grid Search', 2, None, 'grid-search'), + ('Randomized Grid Search', 2, None, 'randomized-grid-search'), ('Optimization, the central part of any Machine Learning ' 'algortithm', 2, @@ -218,56 +219,57 @@ MathJax.Hub.Config({
  • Thursday September 30
  • Searching for Optimal Regularization Parameters \( \lambda \)
  • Grid Search
    • -
  • Optimization, the central part of any Machine Learning algortithm
    • -
  • Revisiting our Logistic Regression case
    • -
  • The equations to solve
    • -
  • Solving using Newton-Raphson's method
    • -
  • Brief reminder on Newton-Raphson's method
    • -
  • The equations
    • -
  • Simple geometric interpretation
    • -
  • Extending to more than one variable
    • -
  • Steepest descent
    • -
  • More on Steepest descent
    • -
  • The ideal
    • -
  • The sensitiveness of the gradient descent
    • -
  • Convex functions
    • -
  • Convex function
    • -
  • Conditions on convex functions
    • -
  • More on convex functions
    • -
  • Some simple problems
    • -
  • Standard steepest descent
    • -
  • Gradient method
    • -
  • Steepest descent method
    • -
  • Steepest descent method
    • -
  • Final expressions
    • -
  • Steepest descent example
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method and iterations
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method
    • -
  • Conjugate gradient method
    • -
  • Revisiting our first homework
    • -
  • Gradient descent example
    • -
  • The derivative of the cost/loss function
    • -
  • The Hessian matrix
    • -
  • Simple program
    • -
  • Gradient Descent Example
    • -
  • And a corresponding example using scikit-learn
    • -
  • Gradient descent and Ridge
    • -
  • Program example for gradient descent with Ridge Regression
    • -
  • Using gradient descent methods, limitations
    • -
  • Friday October 1
    • -
  • Stochastic Gradient Descent
    • -
  • Computation of gradients
    • -
  • SGD example
    • -
  • The gradient step
    • -
  • Simple example code
    • -
  • When do we stop?
    • -
  • Slightly different approach
    • -
  • Program for stochastic gradient
    • +
  • Randomized Grid Search
    • +
  • Optimization, the central part of any Machine Learning algortithm
    • +
  • Revisiting our Logistic Regression case
    • +
  • The equations to solve
    • +
  • Solving using Newton-Raphson's method
    • +
  • Brief reminder on Newton-Raphson's method
    • +
  • The equations
    • +
  • Simple geometric interpretation
    • +
  • Extending to more than one variable
    • +
  • Steepest descent
    • +
  • More on Steepest descent
    • +
  • The ideal
    • +
  • The sensitiveness of the gradient descent
    • +
  • Convex functions
    • +
  • Convex function
    • +
  • Conditions on convex functions
    • +
  • More on convex functions
    • +
  • Some simple problems
    • +
  • Standard steepest descent
    • +
  • Gradient method
    • +
  • Steepest descent method
    • +
  • Steepest descent method
    • +
  • Final expressions
    • +
  • Steepest descent example
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method and iterations
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method
    • +
  • Conjugate gradient method
    • +
  • Revisiting our first homework
    • +
  • Gradient descent example
    • +
  • The derivative of the cost/loss function
    • +
  • The Hessian matrix
    • +
  • Simple program
    • +
  • Gradient Descent Example
    • +
  • And a corresponding example using scikit-learn
    • +
  • Gradient descent and Ridge
    • +
  • Program example for gradient descent with Ridge Regression
    • +
  • Using gradient descent methods, limitations
    • +
  • Friday October 1
    • +
  • Stochastic Gradient Descent
    • +
  • Computation of gradients
    • +
  • SGD example
    • +
  • The gradient step
    • +
  • Simple example code
    • +
  • When do we stop?
    • +
  • Slightly different approach
    • +
  • Program for stochastic gradient
    • @@ -302,7 +304,7 @@ MathJax.Hub.Config({
    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

    -

    Sep 27, 2021

    +

    Sep 28, 2021


    @@ -326,7 +328,7 @@ MathJax.Hub.Config({

  • 9
  • 10
  • ...
    • -
  • 55
    • +
  • 56
  • »
    • diff --git a/doc/pub/week39/html/week39-reveal.html b/doc/pub/week39/html/week39-reveal.html index 75d3f3b38..a7c3d44e7 100644 --- a/doc/pub/week39/html/week39-reveal.html +++ b/doc/pub/week39/html/week39-reveal.html @@ -148,7 +148,7 @@ MathJax.Hub.Config({
    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

     
    -

    Sep 27, 2021

    +

    Sep 28, 2021


    @@ -162,7 +162,7 @@ MathJax.Hub.Config({

    Plan for week 39

    @@ -196,44 +196,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)). 
    import numpy as np
+import pandas as pd
 import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
 
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(315)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
 nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
+MSERidgePredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 2, nlambdas)
+for i in range(nlambdas):
+    lmb = lambdas[i]
+    RegRidge = linear_model.Ridge(lmb)
+    RegRidge.fit(X_train,y_train)
+    ypredictRidge = RegRidge.predict(X_test)
+    MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+# Now plot the results
 plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
 plt.xlabel('log10(lambda)')
 plt.ylabel('MSE')
 plt.legend()
 plt.show()

    -We see from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \). -In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area. +Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function RidgeCV instead of just calling the Ridge function. For RidgeCV we need to passe the array of \( \lambda \) values. +By inspecting the figure we can in turn determine which is the optimal regularization parameter. +This becomes however less functional in the long run. @@ -249,40 +259,109 @@ example here. 

    import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
+from sklearn.model_selection import train_test_split
 from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import GridSearchCV
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
+nlambdas = 10
+lambdas = np.logspace(-4, 2, nlambdas)
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")
+
    +

    +By default the grid search function includes cross validation with five folds. The Scikit-Learn documentation contains more information on how to set the different parameters. + + + +

    + + +

    +An alternative to the above manual grid set up, is to use a random +search where the parameters are tuned from a random distribution +(uniform below) for a fixed number of iterations. A model is +constructed and evaluated for each combination of chosen parameters. +We repeat the previous example but now with a random search. + +

    + + +

    import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import GridSearchCV
+from scipy.stats import uniform as randuniform
+from sklearn.model_selection import RandomizedSearchCV
+
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+param_grid = {'alpha': randuniform()}
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")
    diff --git a/doc/pub/week39/html/week39-solarized.html b/doc/pub/week39/html/week39-solarized.html index 8fdd0accc..a7faaa359 100644 --- a/doc/pub/week39/html/week39-solarized.html +++ b/doc/pub/week39/html/week39-solarized.html @@ -68,6 +68,7 @@ div { text-align: justify; text-justify: inter-word; } None, 'searching-for-optimal-regularization-parameters-lambda'), ('Grid Search', 2, None, 'grid-search'), + ('Randomized Grid Search', 2, None, 'randomized-grid-search'), ('Optimization, the central part of any Machine Learning ' 'algortithm', 2, @@ -238,7 +239,7 @@ MathJax.Hub.Config({
    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

    -

    Sep 27, 2021

    +

    Sep 28, 2021












    @@ -246,7 +247,7 @@ MathJax.Hub.Config({

    Plan for week 39

    @@ -279,44 +280,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)). 
    import numpy as np
+import pandas as pd
 import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
 
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(315)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
 nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
+MSERidgePredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 2, nlambdas)
+for i in range(nlambdas):
+    lmb = lambdas[i]
+    RegRidge = linear_model.Ridge(lmb)
+    RegRidge.fit(X_train,y_train)
+    ypredictRidge = RegRidge.predict(X_test)
+    MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+# Now plot the results
 plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
 plt.xlabel('log10(lambda)')
 plt.ylabel('MSE')
 plt.legend()
 plt.show()

    -We see from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \). -In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area. +Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function RidgeCV instead of just calling the Ridge function. For RidgeCV we need to passe the array of \( \lambda \) values. +By inspecting the figure we can in turn determine which is the optimal regularization parameter. +This becomes however less functional in the long run.











    @@ -332,40 +343,109 @@ example here. 

    import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
+from sklearn.model_selection import train_test_split
 from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import GridSearchCV
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
+nlambdas = 10
+lambdas = np.logspace(-4, 2, nlambdas)
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")
+
    +

    +By default the grid search function includes cross validation with five folds. The Scikit-Learn documentation contains more information on how to set the different parameters. + +

    +









    + +

    + +

    +An alternative to the above manual grid set up, is to use a random +search where the parameters are tuned from a random distribution +(uniform below) for a fixed number of iterations. A model is +constructed and evaluated for each combination of chosen parameters. +We repeat the previous example but now with a random search. + +

    + + +

    import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import GridSearchCV
+from scipy.stats import uniform as randuniform
+from sklearn.model_selection import RandomizedSearchCV
+
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+param_grid = {'alpha': randuniform()}
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")











    diff --git a/doc/pub/week39/html/week39.html b/doc/pub/week39/html/week39.html index 0b5b084ca..6c0ca6be3 100644 --- a/doc/pub/week39/html/week39.html +++ b/doc/pub/week39/html/week39.html @@ -73,6 +73,7 @@ div { text-align: justify; text-justify: inter-word; } None, 'searching-for-optimal-regularization-parameters-lambda'), ('Grid Search', 2, None, 'grid-search'), + ('Randomized Grid Search', 2, None, 'randomized-grid-search'), ('Optimization, the central part of any Machine Learning ' 'algortithm', 2, @@ -243,7 +244,7 @@ MathJax.Hub.Config({

    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

    -

    Sep 27, 2021

    +

    Sep 28, 2021












    @@ -251,7 +252,7 @@ MathJax.Hub.Config({

    Plan for week 39

    @@ -284,44 +285,54 @@ which results in optimal scores (for example the smallest MSE or an \( R2=1 \)). 
    import numpy as np
+import pandas as pd
 import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
-from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import train_test_split
+from sklearn import linear_model
 
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(315)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
 nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
+MSERidgePredict = np.zeros(nlambdas)
+lambdas = np.logspace(-4, 2, nlambdas)
+for i in range(nlambdas):
+    lmb = lambdas[i]
+    RegRidge = linear_model.Ridge(lmb)
+    RegRidge.fit(X_train,y_train)
+    ypredictRidge = RegRidge.predict(X_test)
+    MSERidgePredict[i] = MSE(y_test,ypredictRidge)
+
+# Now plot the results
 plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
+plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
 plt.xlabel('log10(lambda)')
 plt.ylabel('MSE')
 plt.legend()
 plt.show()

    -We see from this plot that the optimal MSE occurs for a value of \( \lambda\in [10,100] \). -In order to nail down the best value of \( \lambda \), we could in turn narrow down the search area. +Here we have performed a rather data greedy calculation as function of the regularization parameter \( \lambda \). There is no resampling here. The latter can easily be added by employing the function RidgeCV instead of just calling the Ridge function. For RidgeCV we need to passe the array of \( \lambda \) values. +By inspecting the figure we can in turn determine which is the optimal regularization parameter. +This becomes however less functional in the long run.











    @@ -337,40 +348,109 @@ example here. 

    import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.model_selection import KFold
+from sklearn.model_selection import train_test_split
 from sklearn.linear_model import Ridge
-from sklearn.model_selection import cross_val_score
-from sklearn.preprocessing import PolynomialFeatures
+from sklearn.model_selection import GridSearchCV
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
 
 # A seed just to ensure that the random numbers are the same for every run.
-np.random.seed(3155)
-# Generate the data.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
 n = 100
-x = np.linspace(-3, 3, n).reshape(-1, 1)
-y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)
-# Decide degree on polynomial to fit
-poly = PolynomialFeatures(degree = 10)
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
 
 # Decide which values of lambda to use
-nlambdas = 500
-lambdas = np.logspace(-3, 5, nlambdas)
-# Initialize a KFold instance
-k = 5
-kfold = KFold(n_splits = k)
-estimated_mse_sklearn = np.zeros(nlambdas)
-i = 0
-for lmb in lambdas:
-    ridge = Ridge(alpha = lmb)
-    estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)
-    estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)
-    i += 1
-plt.figure()
-plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')
-plt.xlabel('log10(lambda)')
-plt.ylabel('MSE')
-plt.legend()
-plt.show()
+nlambdas = 10
+lambdas = np.logspace(-4, 2, nlambdas)
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")
+
    +

    +By default the grid search function includes cross validation with five folds. The Scikit-Learn documentation contains more information on how to set the different parameters. + +

    +









    + +

    + +

    +An alternative to the above manual grid set up, is to use a random +search where the parameters are tuned from a random distribution +(uniform below) for a fixed number of iterations. A model is +constructed and evaluated for each combination of chosen parameters. +We repeat the previous example but now with a random search. + +

    + + +

    import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import GridSearchCV
+from scipy.stats import uniform as randuniform
+from sklearn.model_selection import RandomizedSearchCV
+
+
+def R2(y_data, y_model):
+    return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
+
+def MSE(y_data,y_model):
+    n = np.size(y_model)
+    return np.sum((y_data-y_model)**2)/n
+
+# A seed just to ensure that the random numbers are the same for every run.
+# Useful for eventual debugging.
+np.random.seed(2021)
+
+n = 100
+x = np.random.rand(n)
+y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
+
+Maxpolydegree = 5
+X = np.zeros((n,Maxpolydegree-1))
+
+for degree in range(1,Maxpolydegree): #No intercept column
+    X[:,degree-1] = x**(degree)
+
+# We split the data in test and training data
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
+
+param_grid = {'alpha': randuniform()}
+# create and fit a ridge regression model, testing each alpha
+model = Ridge()
+gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)
+gridsearch.fit(X_train, y_train)
+print(gridsearch)
+ypredictRidge = gridsearch.predict(X_test)
+# summarize the results of the grid search
+print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}")
+print(f"MSE score: {MSE(y_test,ypredictRidge)}")
+print(f"R2 score: {R2(y_test,ypredictRidge)}")











    diff --git a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz index 9b8edb6ef..54492e2af 100644 Binary files a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz and b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz differ diff --git a/doc/pub/week39/ipynb/week39.ipynb b/doc/pub/week39/ipynb/week39.ipynb index 55531bbb4..ffd705dee 100644 --- a/doc/pub/week39/ipynb/week39.ipynb +++ b/doc/pub/week39/ipynb/week39.ipynb @@ -10,7 +10,7 @@ " \n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n", "\n", - "Date: **Sep 27, 2021**\n", + "Date: **Sep 28, 2021**\n", "\n", "Copyright 1999-2021, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n", "\n", @@ -18,7 +18,7 @@ "\n", "## Plan for week 39\n", "\n", - "* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1\n", + "* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression\n", "\n", "* Friday: Stochastic Gradient descent with examples and automatic differentiation\n", "\n", @@ -55,36 +55,45 @@ "%matplotlib inline\n", "\n", "import numpy as np\n", + "import pandas as pd\n", "import matplotlib.pyplot as plt\n", - "from sklearn.model_selection import KFold\n", - "from sklearn.linear_model import Ridge\n", - "from sklearn.model_selection import cross_val_score\n", - "from sklearn.preprocessing import PolynomialFeatures\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn import linear_model\n", "\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", "# A seed just to ensure that the random numbers are the same for every run.\n", - "np.random.seed(3155)\n", - "# Generate the data.\n", + "# Useful for eventual debugging.\n", + "np.random.seed(315)\n", + "\n", "n = 100\n", - "x = np.linspace(-3, 3, n).reshape(-1, 1)\n", - "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n", - "# Decide degree on polynomial to fit\n", - "poly = PolynomialFeatures(degree = 10)\n", + "x = np.random.rand(n)\n", + "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n", + "\n", + "Maxpolydegree = 5\n", + "X = np.zeros((n,Maxpolydegree-1))\n", + "\n", + "for degree in range(1,Maxpolydegree): #No intercept column\n", + " X[:,degree-1] = x**(degree)\n", + "\n", + "# We split the data in test and training data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", "\n", "# Decide which values of lambda to use\n", "nlambdas = 500\n", - "lambdas = np.logspace(-3, 5, nlambdas)\n", - "# Initialize a KFold instance\n", - "k = 5\n", - "kfold = KFold(n_splits = k)\n", - "estimated_mse_sklearn = np.zeros(nlambdas)\n", - "i = 0\n", - "for lmb in lambdas:\n", - " ridge = Ridge(alpha = lmb)\n", - " estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)\n", - " estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n", - " i += 1\n", + "MSERidgePredict = np.zeros(nlambdas)\n", + "lambdas = np.logspace(-4, 2, nlambdas)\n", + "for i in range(nlambdas):\n", + " lmb = lambdas[i]\n", + " RegRidge = linear_model.Ridge(lmb)\n", + " RegRidge.fit(X_train,y_train)\n", + " ypredictRidge = RegRidge.predict(X_test)\n", + " MSERidgePredict[i] = MSE(y_test,ypredictRidge)\n", + "\n", + "# Now plot the results\n", "plt.figure()\n", - "plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n", + "plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')\n", "plt.xlabel('log10(lambda)')\n", "plt.ylabel('MSE')\n", "plt.legend()\n", @@ -95,8 +104,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We see from this plot that the optimal MSE occurs for a value of $\\lambda\\in [10,100]$.\n", - "In order to nail down the best value of $\\lambda$, we could in turn narrow down the search area.\n", + "Here we have performed a rather data greedy calculation as function of the regularization parameter $\\lambda$. There is no resampling here. The latter can easily be added by employing the function **RidgeCV** instead of just calling the **Ridge** function. For **RidgeCV** we need to passe the array of $\\lambda$ values.\n", + "By inspecting the figure we can in turn determine which is the optimal regularization parameter.\n", + "This becomes however less functional in the long run. \n", + "\n", "\n", "## Grid Search\n", "\n", @@ -116,40 +127,116 @@ "outputs": [], "source": [ "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from sklearn.model_selection import KFold\n", + "from sklearn.model_selection import train_test_split\n", "from sklearn.linear_model import Ridge\n", - "from sklearn.model_selection import cross_val_score\n", - "from sklearn.preprocessing import PolynomialFeatures\n", + "from sklearn.model_selection import GridSearchCV\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", "\n", "# A seed just to ensure that the random numbers are the same for every run.\n", - "np.random.seed(3155)\n", - "# Generate the data.\n", + "# Useful for eventual debugging.\n", + "np.random.seed(2021)\n", + "\n", "n = 100\n", - "x = np.linspace(-3, 3, n).reshape(-1, 1)\n", - "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)+ np.random.normal(0, 0.1, x.shape)\n", - "# Decide degree on polynomial to fit\n", - "poly = PolynomialFeatures(degree = 10)\n", + "x = np.random.rand(n)\n", + "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n", + "\n", + "Maxpolydegree = 5\n", + "X = np.zeros((n,Maxpolydegree-1))\n", + "\n", + "for degree in range(1,Maxpolydegree): #No intercept column\n", + " X[:,degree-1] = x**(degree)\n", + "\n", + "# We split the data in test and training data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", "\n", "# Decide which values of lambda to use\n", - "nlambdas = 500\n", - "lambdas = np.logspace(-3, 5, nlambdas)\n", - "# Initialize a KFold instance\n", - "k = 5\n", - "kfold = KFold(n_splits = k)\n", - "estimated_mse_sklearn = np.zeros(nlambdas)\n", - "i = 0\n", - "for lmb in lambdas:\n", - " ridge = Ridge(alpha = lmb)\n", - " estimated_mse_folds = cross_val_score(ridge, x, y, scoring='neg_mean_squared_error', cv=kfold)\n", - " estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n", - " i += 1\n", - "plt.figure()\n", - "plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n", - "plt.xlabel('log10(lambda)')\n", - "plt.ylabel('MSE')\n", - "plt.legend()\n", - "plt.show()" + "nlambdas = 10\n", + "lambdas = np.logspace(-4, 2, nlambdas)\n", + "# create and fit a ridge regression model, testing each alpha\n", + "model = Ridge()\n", + "gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas))\n", + "gridsearch.fit(X_train, y_train)\n", + "print(gridsearch)\n", + "ypredictRidge = gridsearch.predict(X_test)\n", + "# summarize the results of the grid search\n", + "print(f\"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}\")\n", + "print(f\"MSE score: {MSE(y_test,ypredictRidge)}\")\n", + "print(f\"R2 score: {R2(y_test,ypredictRidge)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default the grid search function includes cross validation with five folds. The [Scikit-Learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV) contains more information on how to set the different parameters.\n", + "\n", + "## Randomized Grid Search\n", + "\n", + "An alternative to the above manual grid set up, is to use a random\n", + "search where the parameters are tuned from a random distribution\n", + "(uniform below) for a fixed number of iterations. A model is\n", + "constructed and evaluated for each combination of chosen parameters.\n", + "We repeat the previous example but now with a random search." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.linear_model import Ridge\n", + "from sklearn.model_selection import GridSearchCV\n", + "from scipy.stats import uniform as randuniform\n", + "from sklearn.model_selection import RandomizedSearchCV\n", + "\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "# A seed just to ensure that the random numbers are the same for every run.\n", + "# Useful for eventual debugging.\n", + "np.random.seed(2021)\n", + "\n", + "n = 100\n", + "x = np.random.rand(n)\n", + "y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)\n", + "\n", + "Maxpolydegree = 5\n", + "X = np.zeros((n,Maxpolydegree-1))\n", + "\n", + "for degree in range(1,Maxpolydegree): #No intercept column\n", + " X[:,degree-1] = x**(degree)\n", + "\n", + "# We split the data in test and training data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "\n", + "param_grid = {'alpha': randuniform()}\n", + "# create and fit a ridge regression model, testing each alpha\n", + "model = Ridge()\n", + "gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100)\n", + "gridsearch.fit(X_train, y_train)\n", + "print(gridsearch)\n", + "ypredictRidge = gridsearch.predict(X_test)\n", + "# summarize the results of the grid search\n", + "print(f\"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}\")\n", + "print(f\"MSE score: {MSE(y_test,ypredictRidge)}\")\n", + "print(f\"R2 score: {R2(y_test,ypredictRidge)}\")" ] }, { diff --git a/doc/src/week39/test1.py b/doc/src/week39/codes/test1.py similarity index 100% rename from doc/src/week39/test1.py rename to doc/src/week39/codes/test1.py diff --git a/doc/src/week39/test2.py b/doc/src/week39/codes/test2.py similarity index 100% rename from doc/src/week39/test2.py rename to doc/src/week39/codes/test2.py diff --git a/doc/src/week39/test3.py b/doc/src/week39/codes/test3.py similarity index 100% rename from doc/src/week39/test3.py rename to doc/src/week39/codes/test3.py diff --git a/doc/src/week39/test4.py b/doc/src/week39/codes/test4.py similarity index 100% rename from doc/src/week39/test4.py rename to doc/src/week39/codes/test4.py diff --git a/doc/src/week39/codes/test5.py b/doc/src/week39/codes/test5.py new file mode 100644 index 000000000..7e5830c0c --- /dev/null +++ b/doc/src/week39/codes/test5.py @@ -0,0 +1,43 @@ +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.linear_model import Ridge +from sklearn.model_selection import GridSearchCV + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) + +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(2021) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + +# Decide which values of lambda to use +nlambdas = 10 +lambdas = np.logspace(-4, 2, nlambdas) +# create and fit a ridge regression model, testing each alpha +model = Ridge() +gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas)) +gridsearch.fit(X_train, y_train) +print(gridsearch) +ypredictRidge = gridsearch.predict(X_test) +# summarize the results of the grid search +#print(gridsearch.best_score_) +print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}") +print(f"MSE score: {MSE(y_test,ypredictRidge)}") +print(f"R2 score: {R2(y_test,ypredictRidge)}") diff --git a/doc/src/week39/test6.py b/doc/src/week39/codes/test6.py similarity index 100% rename from doc/src/week39/test6.py rename to doc/src/week39/codes/test6.py diff --git a/doc/src/week39/codes/test7.py b/doc/src/week39/codes/test7.py new file mode 100644 index 000000000..2f8389557 --- /dev/null +++ b/doc/src/week39/codes/test7.py @@ -0,0 +1,44 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn import linear_model + +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(315) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + +# Decide which values of lambda to use +nlambdas = 500 +MSERidgePredict = np.zeros(nlambdas) +lambdas = np.logspace(-4, 2, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + RegRidge = linear_model.Ridge(lmb) + RegRidge.fit(X_train,y_train) + ypredictRidge = RegRidge.predict(X_test) + MSERidgePredict[i] = MSE(y_test,ypredictRidge) + +# Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() diff --git a/doc/src/week39/test5.py b/doc/src/week39/test5.py deleted file mode 100644 index 430ab655c..000000000 --- a/doc/src/week39/test5.py +++ /dev/null @@ -1,61 +0,0 @@ -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from sklearn.model_selection import train_test_split -from sklearn.linear_model import Ridge -from sklearn.preprocessing import StandardScaler -from sklearn.model_selection import cross_val_score -from sklearn.model_selection import GridSearchCV - - - -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -# A seed just to ensure that the random numbers are the same for every run. -# Useful for eventual debugging. -np.random.seed(315) - -n = 100 -x = np.random.rand(n) -y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) - -Maxpolydegree = 5 -X = np.zeros((n,Maxpolydegree-1)) - -for degree in range(1,Maxpolydegree): #No intercept column - X[:,degree-1] = x**(degree) - -# We split the data in test and training data -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - -#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable -X_train_mean = np.mean(X_train,axis=0) -#Center by removing mean from each feature -X_train_scaled = X_train - X_train_mean -X_test_scaled = X_test - X_train_mean -#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered) -#Remove the intercept from the training data. -y_scaler = np.mean(y_train) -y_train_scaled = y_train - y_scaler - -p = Maxpolydegree-1 -I = np.eye(p,p) -# Decide which values of lambda to use -nlambdas = 10 -MSEOwnRidgePredict = np.zeros(nlambdas) -MSERidgePredict = np.zeros(nlambdas) - -lambdas = np.logspace(-4, 2, nlambdas) - -# create and fit a ridge regression model, testing each alpha -model = Ridge() -grid = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas)) -grid.fit(X_train_scaled, y_train_scaled) -print(grid) -ypredictRidge = grid.predict(X_test_scaled) -# summarize the results of the grid search -print(grid.best_score_) -print(grid.best_estimator_.alpha) -print(MSE(y_test,ypredictRidge)) diff --git a/doc/src/week39/week39.do.txt b/doc/src/week39/week39.do.txt index 25f3525dc..816605db6 100644 --- a/doc/src/week39/week39.do.txt +++ b/doc/src/week39/week39.do.txt @@ -5,7 +5,7 @@ DATE: today !split ===== Plan for week 39 ===== -* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 +* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Discussion of project 1 and examples on how to implement Logistic Regression * Friday: Stochastic Gradient descent with examples and automatic differentiation @@ -31,10 +31,57 @@ in evaluating the MSE as function of different $\lambda$ values. Based on these calculations, one tries then to determine the value of the hyperparameter $\lambda$ which results in optimal scores (for example the smallest MSE or an $R2=1$). !bc pycod +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split +from sklearn import linear_model + +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(315) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + +# Decide which values of lambda to use +nlambdas = 500 +MSERidgePredict = np.zeros(nlambdas) +lambdas = np.logspace(-4, 2, nlambdas) +for i in range(nlambdas): + lmb = lambdas[i] + RegRidge = linear_model.Ridge(lmb) + RegRidge.fit(X_train,y_train) + ypredictRidge = RegRidge.predict(X_test) + MSERidgePredict[i] = MSE(y_test,ypredictRidge) + +# Now plot the results +plt.figure() +plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test') +plt.xlabel('log10(lambda)') +plt.ylabel('MSE') +plt.legend() +plt.show() + !ec -We see from this plot that the optimal MSE occurs for a value of $\lambda\in [10,100]$. -In order to nail down the best value of $\lambda$, we could in turn narrow down the search area. +Here we have performed a rather data greedy calculation as function of the regularization parameter $\lambda$. There is no resampling here. The latter can easily be added by employing the function _RidgeCV_ instead of just calling the _Ridge_ function. For _RidgeCV_ we need to passe the array of $\lambda$ values. +By inspecting the figure we can in turn determine which is the optimal regularization parameter. +This becomes however less functional in the long run. + !split ===== Grid Search ===== @@ -45,12 +92,111 @@ included with the library _Scikit-Learn_, as demonstrated for the same example here. !bc pycod +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.linear_model import Ridge +from sklearn.model_selection import GridSearchCV + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) + +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(2021) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + +# Decide which values of lambda to use +nlambdas = 10 +lambdas = np.logspace(-4, 2, nlambdas) +# create and fit a ridge regression model, testing each alpha +model = Ridge() +gridsearch = GridSearchCV(estimator=model, param_grid=dict(alpha=lambdas)) +gridsearch.fit(X_train, y_train) +print(gridsearch) +ypredictRidge = gridsearch.predict(X_test) +# summarize the results of the grid search +print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}") +print(f"MSE score: {MSE(y_test,ypredictRidge)}") +print(f"R2 score: {R2(y_test,ypredictRidge)}") !ec +By default the grid search function includes cross validation with five folds. The "Scikit-Learn documentation":"https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV" contains more information on how to set the different parameters. + +!split +===== Randomized Grid Search ===== + +An alternative to the above manual grid set up, is to use a random +search where the parameters are tuned from a random distribution +(uniform below) for a fixed number of iterations. A model is +constructed and evaluated for each combination of chosen parameters. +We repeat the previous example but now with a random search. +!bc pycod +import numpy as np +from sklearn.model_selection import train_test_split +from sklearn.linear_model import Ridge +from sklearn.model_selection import GridSearchCV +from scipy.stats import uniform as randuniform +from sklearn.model_selection import RandomizedSearchCV + + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) + +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(2021) + +n = 100 +x = np.random.rand(n) +y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + +Maxpolydegree = 5 +X = np.zeros((n,Maxpolydegree-1)) + +for degree in range(1,Maxpolydegree): #No intercept column + X[:,degree-1] = x**(degree) + +# We split the data in test and training data +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + +param_grid = {'alpha': randuniform()} +# create and fit a ridge regression model, testing each alpha +model = Ridge() +gridsearch = RandomizedSearchCV(estimator=model, param_distributions=param_grid, n_iter=100) +gridsearch.fit(X_train, y_train) +print(gridsearch) +ypredictRidge = gridsearch.predict(X_test) +# summarize the results of the grid search +print(f"Best estimated lambda-value: {gridsearch.best_estimator_.alpha}") +print(f"MSE score: {MSE(y_test,ypredictRidge)}") +print(f"R2 score: {R2(y_test,ypredictRidge)}") + +!ec + !split