From c1cae8a2f51665e71e4b50e5f84dcff5d5fe2b15 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Thu, 22 Sep 2022 11:40:55 +0200 Subject: [PATCH] 03: Do the first ones --- Exercise sheet 3/exercise_sheet_03.ipynb | 64 ++++++++++++++++++------ 1 file changed, 49 insertions(+), 15 deletions(-) diff --git a/Exercise sheet 3/exercise_sheet_03.ipynb b/Exercise sheet 3/exercise_sheet_03.ipynb index de0387c..806e1c2 100644 --- a/Exercise sheet 3/exercise_sheet_03.ipynb +++ b/Exercise sheet 3/exercise_sheet_03.ipynb @@ -20,12 +20,12 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "d8b13e80", "metadata": {}, "outputs": [], "source": [ - "NAME = \"\"\n", + "NAME = \"Kees van Kempen\"\n", "NAMES_OF_COLLABORATORS = \"\"" ] }, @@ -62,7 +62,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "948639fc", "metadata": { "deletable": false, @@ -135,7 +135,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "564c7f00", "metadata": { "deletable": false, @@ -154,8 +154,10 @@ "source": [ "def f_inverse_Z(u):\n", " '''Compute the inverse CDF of Z, i.e. F_Z^{-1}(u) for 0 <= u <= 1.'''\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()\n", + " assert 0 <= u and u <= 1\n", + " \n", + " # we need to round to return a natural number (excluding zero, yes).\n", + " return np.ceil(u/(1 - u)).astype(int)\n", "\n", "def random_Z():\n", " return int(f_inverse_Z(rng.random())) # make sure to return an integer" @@ -163,7 +165,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "id": "2ca9d9ca", "metadata": { "deletable": false, @@ -211,7 +213,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "8e201507", "metadata": { "deletable": false, @@ -226,24 +228,56 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdY0lEQVR4nO3dfZzNdf7/8cfLGKELylWicpEbUYhBKEPC1ipMSt1ascWkJL7f32qlzVarbb9xayuFrK1QbRYp0m6bXbNqU7loulAhF635KmQzMYZmeP3+mJnzNcyM+ejMfM6ced5vt7lxznwuXuczc+Z5Xp+L98fcHRERkdKqEnYBIiJSsSg4REQkEAWHiIgEouAQEZFAFBwiIhJI1bALKA9169b1Jk2ahF2GiEiFsXbt2m/dvV5R34vr4DCza4BrLrjgAtasWRN2OSIiFYaZfVXc9+J6V5W7L3X31Fq1aoVdiohI3Ijr4BARkehTcIiISCBxfYxDyl9OTg4ZGRkcPHgw7FJEpBSqV69O48aNSUxMLPU8Cg6JqoyMDE4//XSaNGmCmYVdjoiUwN3Zs2cPGRkZNG3atNTzxfWuKjO7xsxmZWZmhl1KpXHw4EHq1Kmj0BCpAMyMOnXqBN5DENfBobOqwqHQEKk4Tub9GtfBISIi0adjHCfwyoavwy4BgJSWDcMu4aREe/sF3Q4PPPAAp512Gr/4xS+iVsPVV1/NSy+9BMBLL73EnXfeGZXljh8/njfeeIOrr76aKVOmRGWZP8Zvf/tbJk6cGHYZJVqzZg1z587lySefDDzv448/TmpqKjVr1gT+7+dau3btKFdZdtLS0pg6dSqvv/56ua5XHYdIQG+88Qa1a9dm7969TJ8+PWrLfeaZZ1i3bl1MhAbkBUcsyM3NLfZ7SUlJJxUakBccBw4ciDwu+LnKianjCGDS0OvKdX0PzVtUruuLFw8//DBz587l3HPPpV69enTs2BGAzZs3M3r0aHbv3k3NmjX5wx/+QKtWrRg+fDhnnHEGa9as4ZtvvuHRRx9l8ODBfP311wwZMoTvv/+e3NxcZsyYweWXX06TJk1Ys2YNEyZMYPPmzbRv354+ffrwzTffMHjwYAYMGADAzTffzJAhQ7j22msjtbk799xzD3/5y18wM371q19FpsnKyqJLly7ce++9DBkyJDLPBx98wLhx48jOzqZGjRo899xztGzZksOHD/PLX/6SN998EzNj5MiRjBkzhtWrVzN27FiysrI45ZRT+Pvf/07NmjWZMGECaWlpHDp0iNGjR3P77beTlpbGpEmTqFOnDhs2bKBHjx5Mnz6diRMnkp2dTfv27WnTpg0vvvgiAwcOZPv27Rw8eJCxY8eSmpoKwGmnncbYsWN5/fXXqVGjBq+99hoNGjRg586djBo1ii1btgAwY8YMunXrxgsvvMCTTz7JDz/8QJcuXZg+fToJCQmFfobPP/88y5Yt4+DBg2RlZbF06VLGjBnDJ598Qm5uLg888AADBgwo9Ik7KyuryGmK2k7uzo4dO+jVqxd169ZlxYoVkZ9r3bp1eeyxx3j22WcBGDFiBOPGjWPbtm1cddVVXHbZZbz77rs0atSI1157jRo1ahSqfcGCBTz44IMkJCRQq1YtVq5cybZt2xg6dChZWVkAPPXUU3Tr1o20tDR+/etf06BBA9LT00lJSeHiiy/miSeeIDs7m1dffZXmzZszfPhwqlevzvr169m5cyePPfYY/fv3L7Te4l7/+vXr+fnPf84PP/zAkSNHWLRoES1atPhxbzJ3j/uvjh07+sla9MWOyFebTl3L9evodVcUn332WaHHR7+GaHydyJo1a/yiiy7yrKwsz8zM9ObNm/uUKVPc3f2KK67wjRs3urv7e++957169XJ392HDhvngwYP98OHDvn79em/evLm7u0+dOtUnT57s7u65ubn+/fffu7v7+eef77t37/atW7d6mzZtIutOS0vzAQMGuLv73r17vUmTJp6Tk1OovoULF/qVV17pubm5/s033/i5557rO3bkva5TTz21yNeUmZkZWc5bb73lKSkp7u4+ffp0T0lJiXxvz549fujQIW/atKl/8MEHheZ95pln/De/+Y27ux88eNA7duzoW7Zs8RUrVvgpp5zimzdv9tzcXL/yyit9wYIFRdazZ88ed3c/cOCAt2nTxr/99lt3dwd8yZIl7u4+fvz4yHpuuOEG//3vfx/Zfnv37vXPPvvM+/fv7z/88IO7u99xxx0+Z86c417zc889540aNYqs89577/V58+a5u/t3333nLVq08P379/uKFSv8pz/9aYnTFLWdjv45Fih4XPA7tH//ft+3b5+3bt3a161b51u3bvWEhAT/8MMP3d39+uuvj6zvaBdddJFnZGRE6nB3z8rK8uzsbHd337hxoxf8TVqxYoXXqlXLd+zY4QcPHvRzzjnHJ02a5O7ujz/+uI8dO9bd835H+/Xr54cPH/aNGzd6o0aNPDs7u1Sv/6677vIXXnjB3d0PHTrkBw4cOK7mY9+37u7AGi/mb2pcdxxHD3IYDeoAYt/bb7/NoEGDIvutCz7t79+/n3fffZfrr78+Mu2hQ4ci/x84cCBVqlShdevW7Ny5E4BOnTpx6623kpOTw8CBA2nfvn2J605OTmb06NHs2rWLV155heuuu46qVQu/xd555x1uuukmEhISaNCgAcnJyaxevbpQV3KszMxMhg0bxqZNmzAzcnJyAFi+fDmjRo2KrOOss87ik08+oWHDhnTq1AmAM844A4C//e1vfPzxxyxcuDCyzE2bNlGtWjU6d+5Ms2bNALjpppt45513GDx48HF1PPnkkyxevBiA7du3s2nTJurUqUO1atUin347duzIW2+9BcA//vEP5s6dCxD59D1v3jzWrl0bqS87O5v69esX+br79OnDWWedFal/yZIlTJ06Fcg77fvf//53oemLm6ao7VSSd955h0GDBnHqqacCkJKSwttvv821115L06ZNI78HHTt2ZNu2bcfN3717d4YPH84NN9xASkoKkHdh7F133UV6ejoJCQls3LgxMn2nTp1o2DDv2F3z5s3p27cvABdffDErVqyITHfDDTdQpUoVWrRoQbNmzfjiiy9K9fq7du3Kww8/TEZGBikpKT++2yDOd1W5+1JgaVJS0siwa5HyU9TphUeOHKF27dqkp6cXOc8pp5wS+X/ehy3o0aMHK1euZNmyZQwdOpTx48dzyy23lLjuoUOH8uKLL/Lyyy9HdnUcrWDZQdx///306tWLxYsXs23bNnr27BlZ1rGvtajnCp6fNm0a/fr1K/R8WlracdMXNX9aWhrLly9n1apV1KxZk549e0bO/U9MTIzMk5CQUOIxCXdn2LBhPPLII4WeX7x4MQ8++CAAs2fPBoj84S6Yb9GiRbRs2bLQfAUhX9I0xW2TkmosztG/JwkJCWRnZx83zcyZM3n//fdZtmwZ7du3Jz09nWnTptGgQQM++ugjjhw5QvXq1YtcZpUqVSKPq1SpUmhbnujnVNzrv/DCC+nSpQvLli2jX79+zJ49myuuuKKkTXBCOjgucaVHjx4sXryY7Oxs9u3bx9KlS4G8T95NmzZlwYIFQN6b7KOPPipxWV999RX169dn5MiR3Hbbbaxbt67Q908//XT27dtX6Lnhw4fz+OOPA9CmTZsi65s/fz6HDx9m9+7drFy5ks6dO5dYR2ZmJo0aNQLy9v0X6Nu3LzNnzoz8cfnPf/5Dq1at2LFjB6tXrwZg37595Obm0q9fP2bMmBHpVjZu3BjZ3/7BBx+wdetWjhw5wvz587nsssuAvEAomD4zM5MzzzyTmjVr8sUXX/Dee++VWDNA7969mTFjBgCHDx/m+++/p3fv3ixcuJBdu3ZFav7qq68YNGgQ6enppKenk5SUdNyy+vXrx7Rp0yJ/1D/88MNST1PUdoKif36Q9zN69dVXOXDgAFlZWSxevJjLL7/8hK+3wObNm+nSpQsPPfQQdevWZfv27WRmZtKwYUOqVKnCvHnzOHz4cKmXV2DBggUcOXKEzZs3s2XLluMCorjXv2XLFpo1a8bdd9/Ntddey8cffxx43ceK645DwlfepxF36NCBIUOG0L59e84///xCb/gXX3yRO+64g8mTJ5OTk8ONN95Iu3btil1WWloaU6ZMITExkdNOOy2y26VAnTp16N69OxdddBFXXXUVU6ZMoUGDBlx44YUMHDiwyGUOGjSIVatW0a5dO8yMRx99lLPPPrvE13TPPfcwbNgwHnvssUKfFEeMGMHGjRtp27YtiYmJjBw5krvuuov58+czZsyYyMH05cuXM2LECLZt20aHDh1wd+rVq8err74KQNeuXZkwYQKffPIJPXr0YNCgQQCkpqbStm1bOnTowLPPPsvMmTNp27YtLVu25NJLLy2xZoAnnniC1NRU/vjHP5KQkMCMGTPo2rUrkydPpm/fvhw5coTExESefvppzj///BKXdf/99zNu3Djatm2Lu9OkSZPIKagFn7yLm6a47ZSamspVV11Fw4YNC+0S6tChA8OHD48E+ogRI7jkkkuK3C1VlPHjx7Np0ybcnd69e9OuXTvuvPNOrrvuOhYsWECvXr0KdVOl1bJlS5KTk9m5cyczZ84s1LWU9Prnz5/PCy+8QGJiImeffTaTJk0KvO5j2cm0zhVNUlKSn+yNnHQdRzCff/45F154YdhlhObAgQNcfPHFrFu3joowYkFY1wFEy6JFi1iyZAlz5swJu5QyNXz4cPr371/ksadoKOp9a2Zr3f349g/tqhKJmuXLl9OqVSvGjBlTIUKjoluyZAn33Xcft99+e9ilVDrqOE5AHUcwlb3jEKmI1HGIiEiZiuvg0LDqIiLRF9fB4RpWXUQk6uI6OEREJPoUHBJ3jh61Ni0t7bjB4Mra888/z44dOyKPR4wYwWeffRZ4OWHULlIaCg6JO9Ee7rwoJQ2rcWxwzJ49m9atW5dpPSLlScEhcefo4c7Hjx/P/v37GTx4MK1ateLmm2+ODMmwdu1akpOT6dixI/369ePrr/NOvU5PT+fSSy+lbdu2DBo0iO+++w6Anj17MnHiRJKTk3niiSeKnH/hwoWsWbOGm2++mfbt25OdnU3Pnj0pOB38r3/9Kx06dKBdu3b07t0byBvyo1u3blxyySV069aNDRs2hLDVREpPQ45I3Pnd737Hp59+Snp6OmlpaZF7Epxzzjl0796df/3rX3Tp0oUxY8bw2muvUa9ePebPn899993Hs88+yy233MK0adNITk5m0qRJPPjgg5Hxp/bu3cs///lPcnJySE5OLnL+p556iqlTpx435tLu3bsZOXIkK1eupGnTppExk1q1asXKlSupWrUqy5cvZ+LEiSxapJGYJXYpOCSqhg0bFhnSvDykpaWdcJrOnTvTuHFjANq3b8+2bduoXbs2n376KX369AHyBuFr2LAhmZmZ7N27l+TkZCDv9Rw9FHvBDZY2bNhQ5Pwlee+99+jRowdNmzYF/m947+KGTReJVQoOiXvHDoWdm5uLu9OmTRtWrVpVaNoTXfNTMDhdcfOXpLjhvYsbNl0kVik4JKrmzJkT+pAjxQ2XfbSWLVuye/duVq1aRdeuXcnJyWHjxo20adOGM888k7fffpvLL7+cefPmRbqP0s5f3Pq7du3K6NGj2bp1a2RX1VlnnVXssOkisSqugyPadwCUiuHo4c5r1KhBgwYNjpumWrVqLFy4kLvvvpvMzExyc3MZN24cbdq0Yc6cOYwaNYoDBw7QrFkznnvuuUDzDx8+nFGjRlGjRo1CHUm9evWYNWsWKSkpHDlyhPr16/PWW28VO2y6SKzSIIcnoEEOg9EghyIVjwY5FBGRMqXgEBGRQBQcEnWVYfenSLw4mfergkOiqnr16uzZs0fhIVIBuDt79uw57v7lJxLXZ1VJ+WvcuDEZGRns3r077FJEpBSqV68euUC2tBQcElWJiYmRK6NFJD5pV5WIiASi4BARkUAUHCIiEoiCQ0REAonr4DCza8xs1olGPBURkdKL6+Bw96XunlqrVq2wSxERiRtxHRwiIhJ9Cg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCievg0B0ARUSiL66DQ3cAFBGJvrgODhERiT4Fh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoFUDbsAKZ1XNnwddgkApLRsGHYJIhIydRwiIhKIOo4KaNLQ68p1fQ/NW1Su6xOR2KaOQ0REAlHHUQGpAxCRMKnjEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQmkwgWHmTUzsz+a2cKwaxERqYzKNTjM7Fkz22Vmnx7z/E/MbIOZfWlmE0pahrtvcffbyrZSEREpTnnfj+N54ClgbsETZpYAPA30ATKA1Wa2BEgAHjlm/lvdfVf5lCoiIkUp1+Bw95Vm1uSYpzsDX7r7FgAzexkY4O6PAP1Pdl1mlgqkApx33nknuxgRETlGLBzjaARsP+pxRv5zRTKzOmY2E7jEzO4tbjp3n+XuSe6eVK9evehVKyJSycXCrWOtiOe8uIndfQ8wquzKERGRksRCx5EBnHvU48bAjpBqERGRE4iF4FgNtDCzpmZWDbgRWBJyTSIiUozAwWFmp+afCRWYmf0JWAW0NLMMM7vN3XOBu4A3gc+BP7v7+pNZfhHru8bMZmVmZkZjcSIiQimOcZhZFfK6gJuBTsAh4BQz2w28Acxy902lWZm731TM82/kLyuq3H0psDQpKWlktJctIlJZlabjWAE0B+4Fznb3c929PnA58B7wOzP7WRnWKCIiMaQ0Z1Vd6e45xz7p7v8BFgGLzCwx6pWJiEhMOmHHURAaZjb52O8VHOsoKlhERCQ+BTk43sjMIscozKw+sDz6JUWPDo6LiERfkOC4HUg1s85m1gn4BzC1bMqKDndf6u6ptWrVCrsUEZG4UZqzquYC64APgdHAS0AuMNDdvyzb8kREJNaUpuOYkz/dreSFRhPgO+BnZja47EoTEZFYdMKOw93/Dvy94LGZVQVaA+2ASwHdUKkSeWXD12GXQErLhmGXIFKplWZXlbl7ZNDB/Cu9P87/mlfUNCIiEr9KdQGgmY0xs0I3tTCzamZ2hZnNAYaVTXk/js6qEhGJvtJcAPgT8o5v/MnMmgJ7gerk3aHvb8Dv3T29rAr8MTTkSNmZNPS6cl3fQ/MWlev6RKR4pTnGcRCYDkzPv0K8LpDt7nvLuDYREYlBgW7k5O45ZnaLu/9PWRUkFYM6AJHK62Tux/GlmT2ef4zjZ2amvyAiIpXIyQTHYvI6lR35/+paDhGRSuRkguNPwF+A7kBf4MyoViQiIjEtcHC4+xB3X+buG4D/Ap6OflnRodNxRUSir9TBYWZXmdn7ZrbBzP5sZl3dfScwogzr+1E0yKGISPQF6TimA/9N3jAjs4ApZnaTu2eVSWUiIhKTgpyOu9Pd/5X//+Vmtgp4n7xjHiIiUkkE6Ti2mdlkM6uW/zgH2FcGNYmISAwLEhwOpADbzewd4EsgzcxalEllIiISk0q9q8rdbwIws+rAReQNq94OmG1mzdz93LIpUUREYkmgIUcgMnbVmvwvERGpZE7mAsAKQ9dxiIhEX1wHh67jEBGJvrgODhERiT4Fh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQkk8JXjImF7ZcPXYZcAQErLhmGXIBKKuO44dOW4iEj0xXXH4e5LgaVJSUkjw65FysakodeV6/oemreoXNcnEoviuuMQEZHoi+uOQ+KfOgCR8qeOQ0REAlFwiIhIIAoOEREJRMEhIiKBKDhERCQQBYeIiASi4BARkUAUHCIiEoiCQ0REAonr4NAghyIi0RfXweHuS909tVatWmGXIiISN+I6OEREJPoUHCIiEoiCQ0REAlFwiIhIILofh8hJ0r3PpbJSxyEiIoGo4xCJAt37XCoTdRwiIhKIOg6RKFAHIJWJOg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISSFwHh27kJCISfXEdHLqRk4hI9MV1cIiISPQpOEREJBAFh4iIBKKxqkQqON0XRMqbOg4REQlEHYdIHNF9QaQ8qOMQEZFA1HGIxBF1AFIe1HGIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggunJcRKJCo/RWHuo4REQkEHUcIhJ1GqU3vqnjEBGRQNRxiEjUqQOIb+o4REQkEAWHiIgEouAQEZFAFBwiIhJIhQsOMxtoZn8ws9fMrG/Y9YiIVDblGhxm9qyZ7TKzT495/idmtsHMvjSzCSUtw91fdfeRwHBgSBmWKyIiRSjv03GfB54C5hY8YWYJwNNAHyADWG1mS4AE4JFj5r/V3Xfl//9X+fOJiEg5KtfgcPeVZtbkmKc7A1+6+xYAM3sZGODujwD9j12GmRnwO+Av7r6uuHWZWSqQCnDeeedF5wWISMyLhTGz4n28rFi4ALARsP2oxxlAlxKmHwNcCdQyswvcfWZRE7n7LGAWQFJSkkepVhGpQDT0SdmIheCwIp4r9g+9uz8JPFl25YiISEliITgygHOPetwY2BFSLSISRypLB1DeYuF03NVACzNrambVgBuBJSHXJCIixSjv03H/BKwCWppZhpnd5u65wF3Am8DnwJ/dfX2U1neNmc3KzMyMxuJERITyP6vqpmKefwN4owzWtxRYmpSUNDLayxYRqaxiYVeViIhUIAoOEREJRMEhIiKBxMLpuGXGzK4BrrngggvCLkVEKpFYuHodyu4K9rjuONx9qbun1qpVK+xSRETiRlx3HCIiYYvHYU/iuuMQEZHoU8chIlKG4nHYE3UcIiISSFwHh4YcERGJvrgODp1VJSISfXEdHCIiEn0KDhERCUTBISIigSg4REQkEAWHiIgEYu4edg1lzsx2A1+FXUcMqAt8G3YRMULbojBtj8K0PeB8d69X1DcqRXBIHjNb4+5JYdcRC7QtCtP2KEzbo2TaVSUiIoEoOEREJBAFR+UyK+wCYoi2RWHaHoVpe5RAxzhERCQQdRwiIhKIgkNERAJRcMQ5MzvXzFaY2edmtt7MxoZdUywwswQz+9DMXg+7ljCZWW0zW2hmX+T/jnQNu6Ywmdl/5b9PPjWzP5lZ9bBrikUKjviXC/w/d78QuBQYbWatQ64pFowFPg+7iBjwBPBXd28FtKMSbxMzawTcDSS5+0VAAnBjuFXFJgVHnHP3r919Xf7/95H3h6FRuFWFy8waAz8FZoddS5jM7AygB/BHAHf/wd33hlpU+KoCNcysKlAT2BFyPTFJwVGJmFkT4BLg/ZBLCdvjwD3AkZDrCFszYDfwXP5uu9lmdmrYRYXF3f8XmAr8G/gayHT3v4VbVWxScFQSZnYasAgY5+7fh11PWMysP7DL3deGXUsMqAp0AGa4+yVAFjAh3JLCY2ZnAgOApsA5wKlm9rNwq4pNCo5KwMwSyQuNF939lbDrCVl34Foz2wa8DFxhZi+EW1JoMoAMdy/oQBeSFySV1ZXAVnff7e45wCtAt5BrikkKjjhnZkbePuzP3f2xsOsJm7vf6+6N3b0JeQc+/+HulfJTpbt/A2w3s5b5T/UGPguxpLD9G7jUzGrmv296U4lPFihJ1bALkDLXHRgKfGJm6fnPTXT3N8IrSWLIGOBFM6sGbAF+HnI9oXH3981sIbCOvLMRP0RDjxRJQ46IiEgg2lUlIiKBKDhERCQQBYeIiASi4BARkUAUHCIiEoiCQ0REAlFwiIhIIAoOkRCY2ZVmNi/sOkROhoJDJBztyLsyWaTCUXCIhKMd8KGZnWJmz5vZb/PHRxKJeRqrSiQc7YBdwJvAbHevrCP0SgWksapEyln+MPffAl8Bt7v7qpBLEglEu6pEyl9rYDV5I7AeDrkWkcAUHCLlrx3wLnn3A3nOzBqEXI9IIAoOkfLXDvjU3TcCvwT+nL/7SqRC0DEOEREJRB2HiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigfx/vf8Mt7nTGkQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "def accept_probability_X(k):\n", " '''Return the appropriate acceptance probability on the event Z=k.'''\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()\n", + " \n", + " # The resulting c is found by doing algebraic things, and seeing that\n", + " # 6/pi^2*(1 + 1/k) <= c should hold for all allowed k. For k = 1, c is\n", + " # the largest, so: tada.\n", + " c = 12*np.pi**(-2)\n", + " \n", + " # But we do not need this variable, we only need the acceptance probability.\n", + " return .5*(1 + 1/k)\n", " \n", "def random_X():\n", " return sample_acceptance_rejection(random_Z,accept_probability_X)\n", "\n", "# Verify numerically\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "samples = [random_X() for _ in range(1000000)]\n", + "k = np.array(range(1,11))\n", + "\n", + "# The theoretical densities to compare the densities from the\n", + "# acceptance-rejection densities to:\n", + "p_X = lambda k: 6/(np.pi*k)**2\n", + "\n", + "plt.figure()\n", + "plt.hist(samples, k - .5, density=True, color=\"lightblue\", label=\"density of acceptance-rejection samples\")\n", + "plt.scatter(k[:-1], [p_X(i) for i in k[:-1]], s=800, marker=\"_\", color=\"black\", label=\"theoretical\")\n", + "plt.yscale(\"log\")\n", + "plt.ylabel(\"$p_X(k)$\")\n", + "plt.xlabel(\"$k$\")\n", + "plt.legend()\n", + "plt.show()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "1b0e1d47", "metadata": { "deletable": false, @@ -581,7 +615,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" },