From 5d75b5c59c3efa792255ea971ceade35469e3165 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Mon, 19 Sep 2022 16:29:52 +0200 Subject: [PATCH] 02: Simulate Z_n from the Perato distribution --- Exercise sheet 2/exercise_sheet_02.ipynb | 35 +++++++++++++++++++++--- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/Exercise sheet 2/exercise_sheet_02.ipynb b/Exercise sheet 2/exercise_sheet_02.ipynb index 95bfd86..5a4e28b 100644 --- a/Exercise sheet 2/exercise_sheet_02.ipynb +++ b/Exercise sheet 2/exercise_sheet_02.ipynb @@ -746,7 +746,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "id": "b06896e5", "metadata": { "deletable": false, @@ -762,12 +762,39 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhk0lEQVR4nO3de3CU9dn/8ffFqakpqFXR/jiU2AdLI8StDfE4slal0arU6eOAsfX06Erroc5Uf6U+1plap8XRaa3W58Eda/VpTeE3KpaxVK2HiCcqaPcJB8WmgCUghaJiUCKi1++P3YSbZZPcSXazm3s/rxmGvQ/f3e/e4pVvru/3vm5zd0REJLqGFLsDIiJSWAr0IiIRp0AvIhJxCvQiIhGnQC8iEnHDit2BXA4++GCfMGFCsbshIjJovPLKK/9y90NyHSvJQD9hwgSWL19e7G6IiAwaZvZmV8dCpW7MrN7M1phZi5nN6ea8qWb2sZn9e2/biohIYfQY6M1sKHAXcDpQDZxnZtVdnHcL8Hhv24qISOGEGdHXAS3uvtbddwHzgRk5zrsKeAjY0oe2IiJSIGFy9GOADYHtVuCY4AlmNgY4B/gqMLU3bQPvkQASAOPHjw/RLREpdR999BGtra20t7cXuyuRUVFRwdixYxk+fHjoNmECveXYl10g53bgB+7+sdlep4dpm97pngSSALW1tSrAIxIBra2tjBw5kgkTJpAVG6QP3J1t27bR2tpKVVVV6HZhAn0rMC6wPRbYlHVOLTA/8x/yYOAMM9sdsq2IRFR7e7uCfB6ZGQcddBBbt27tVbswgX4ZMNHMqoCNwCygIXiCu3f+aDGz+4BH3f0RMxvWU1sRiTYF+fzqy/XsMdC7+24zu5L0apqhwL3uvsrMZmeOz+tt2173UkRE+izUDVPuvhhYnLUvZ4B394t6aisiElWpVIpNmzZxxhln9KpdPB7ntttuo7a2Nu99Ksk7Y2XwmDDnj52v18/9ehF70nuDue9SulKpFMuXL+91oC8kFTUTkUj73e9+R11dHbFYjMsvv5y//OUv1NTU0N7ezvvvv8+RRx7JypUraWpq4qSTTuKcc86hurqa2bNn88knnwDwxBNPcNxxx3H00Udz7rnnsmPHDgCWLVvG8ccfz1FHHUVdXR3bt2/nxhtvZMGCBcRiMRYsWMD777/PJZdcwtSpU/nyl7/MH/7wBwB27tzJrFmzqKmpYebMmezcubNg10AjehEZENdccw2pVCqv7xmLxbj99tu7PP7aa6+xYMECXnjhBYYPH853v/td1qxZw9lnn80NN9zAzp07+da3vsXkyZNpamri5ZdfZvXq1Xz+85+nvr6ehx9+mHg8zs0338yTTz5JZWUlt9xyCz//+c+ZM2cOM2fOZMGCBUydOpX33nuP/fbbj5tuuonly5fzq1/9CoDrr7+er371q9x77728++671NXVceqpp3L33Xez33770dzcTHNzM0cffXRer02QAr2IRNZTTz3FK6+8wtSp6fs4d+7cyejRo7nxxhuZOnUqFRUV3HHHHZ3n19XVcfjhhwNw3nnn8fzzz1NRUcHq1as54YQTANi1axfHHXcca9as4XOf+1zne48aNSpnH5544gkWLVrEbbfdBqSXnP7jH/9gyZIlXH311QDU1NRQU1NTmIuAAr2IDJDuRt6F4u5ceOGF/OxnP9tr/+bNm9mxYwcfffQR7e3tVFZWAvsuXTQz3J3TTjuN3//+93sda25uDrXU0d156KGH+OIXv7jPsYFaeqocvYhE1imnnMKDDz7Ili3pElxvv/02b775JolEgp/85Cecf/75/OAHP+g8/+WXX2bdunV88sknLFiwgBNPPJFjjz2WF154gZaWFgA++OAD3njjDSZNmsSmTZtYtmwZAG1tbezevZuRI0fS1tbW+Z5f+9rXuPPOO3FP3/D/17/+FYCTTjqJBx54AICVK1fS3NxcsOugEb2IRFZ1dTU333wz06dP55NPPmH48OHMmDGDYcOG0dDQwMcff8zxxx/P008/zZAhQzjuuOOYM2cOK1as6JyYHTJkCPfddx/nnXceH374IQA333wzRxxxBAsWLOCqq65i586dfPrTn+bJJ5/k5JNPZu7cucRiMX74wx/yox/9iGuuuYaamhrcnQkTJvDoo4/yne98h4svvpiamhpisRh1dXUFuw7W8VOmlNTW1roePDI4DOYlioO574PFa6+9xpe+9KVidyOUpqYmbrvtNh599NFid6VHua6rmb3i7jkX4St1IyIScUrdiIiQvjM1Ho8XuxsFoRG9iEjEKdCLiEScAr2ISMQp0IuIRJwmY0VkwASXtOZDmGWx69ev58wzz2TlypV77b/xxhs56aSTOPXUU3O2e+SRRzjiiCOorq7OS1+LSSN6ESlLN910U5dBHtKBfvXq1Xn5rN27d+flffpKgV5EIu/jjz/msssu48gjj2T69Ons3LmTiy66iAcffBCAOXPmUF1dTU1NDddeey0vvvgiixYt4rrrriMWi/H3v/+dVCrFscceS01NDeeccw7vvPMOkC5VXFNTw3HHHcd1113H5MmTAbjvvvs499xzOeuss5g+fTo7duzglFNO4eijj2bKlCmd5YrXr1/PpEmTuPTSS5k8eTLnn38+Tz75JCeccAITJ07k5Zdf7vf3V6AXkcj729/+xhVXXMGqVas44IADeOihhzqPvf322yxcuJBVq1bR3NzMDTfcwPHHH8/ZZ5/NrbfeSiqV4gtf+AIXXHABt9xyC83NzUyZMoUf//jHAFx88cXMmzePl156iaFDh+71uS+99BL3338/Tz/9NBUVFSxcuJBXX32VZ555hu9///ud9W9aWlr43ve+R3NzM6+//jqNjY08//zz3Hbbbfz0pz/t9/cPFejNrN7M1phZi5nNyXF8hpk1m1nKzJab2YmBY+vNbEXHsX73WESkl6qqqojFYgB85StfYf369Z3HRo0aRUVFBZdeeikPP/ww++233z7tt2/fzrvvvsu0adMAuPDCC1myZAnvvvsubW1tHH/88QA0NDTs1e60007js5/9LJCuYnn99ddTU1PDqaeeysaNG/nnP//Z2b8pU6YwZMgQjjzySE455RTMjClTpuzV177qMdCb2VDgLuB0oBo4z8yyZyeeAo5y9xhwCXBP1vGT3T3WVR0GEZFC+tSnPtX5eujQoXvlzIcNG8bLL7/MN7/5TR555BHq6+tDv29PtcI6yh8DPPDAA2zdupVXXnmFVCrFoYceSnt7+z79GzJkSOf2kCFD8pLfDzOirwNa3H2tu+8C5gMzgie4+w7f840rgdKrlCYiksOOHTvYvn07Z5xxBrfffnvnU7CC5Yb3339/DjzwQJ577jkAfvvb3zJt2jQOPPBARo4cydKlSwGYP39+l5+zfft2Ro8ezfDhw3nmmWd48803C/vFAsIsrxwDbAhstwLHZJ9kZucAPwNGA8E1Tw48YWYO3O3uyVwfYmYJIAEwfvz4UJ0XkcGlFKuEtrW1MWPGDNrb23F3fvGLXwAwa9YsLrvsMu644w4efPBB7r//fmbPns0HH3zA4Ycfzm9+8xsAfv3rX3PZZZdRWVlJPB5n//33z/k5559/PmeddRa1tbXEYjEmTZo0YN+xxzLFZnYu8DV3vzSz/W2gzt2v6uL8k4Ab3f3UzPb/cfdNZjYa+DNwlbsv6e4zVaZ48BjMpX4Hc98Hi8FUprivduzYwWc+8xkA5s6dy1tvvcUvf/nLgn5mIcoUtwLjAttjgU1dnZwJ4l8ws4Mz25syf28BFpJOBYmIRMIf//hHYrEYkydP5rnnnuOGG24odpf2ESZ1swyYaGZVwEZgFrDX1LKZ/Rvwd3d3MzsaGAFsM7NKYIi7t2VeTwduyus3EBEpopkzZzJz5sxid6NbPQZ6d99tZlcCjwNDgXvdfZWZzc4cnwd8E7jAzD4CdgIzM0H/UGBh5gG4w4BGd3+sQN9FREqQuw/YQ7DLQV+eChiq1o27LwYWZ+2bF3h9C3BLjnZrgaN63SspG93VPsl33jzfdVakZxUVFWzbto2DDjpIwT4P3J1t27ZRUVHRq3YqaiYiBTN27FhaW1vZunVrsbsSGRUVFYwdO7ZXbRToRaRghg8fTlVVVbG7UfZU60ZEJOIU6EVEIk6BXkQk4hToRUQiToFeRCTiFOhFRCJOgV5EJOIU6EVEIk43TMmgpBLDIuFpRC8iEnEK9CIiEadALyIScQr0IiIRp0AvIhJxCvQiIhEXKtCbWb2ZrTGzFjObk+P4DDNrNrOUmS03sxPDthURkcLqMdCb2VDgLuB0oBo4z8yqs057CjjK3WPAJcA9vWgrIiIFFGZEXwe0uPtad98FzAdmBE9w9x2+54m1lYCHbSsiIoUVJtCPATYEtlsz+/ZiZueY2evAH0mP6kO3FRGRwgkT6HM9ut332eG+0N0nAd8AftKbtgBmlsjk95frQcIiIvkTJtC3AuMC22OBTV2d7O5LgC+Y2cG9aevuSXevdffaQw45JES3REQkjDCBfhkw0cyqzGwEMAtYFDzBzP7NzCzz+mhgBLAtTFsRESmsHqtXuvtuM7sSeBwYCtzr7qvMbHbm+Dzgm8AFZvYRsBOYmZmczdm2QN9FRERyCFWm2N0XA4uz9s0LvL4FuCVsWxERGTiqRy8lSzXnRfJDJRBERCJOgV5EJOIU6EVEIk6BXkQk4jQZK70SnCDtr2QyyebGO8Ode/hGEolErz8jbH818StRphG9DKhkMkk8Hicej3P55Zfz4YaVPbbZtWUdjY2NA9A7kWjSiF4GVGNjI6lUilgsxrRp01hRMYWRsfpu22xunEMqlSIej3fuawvRTkTSFOhlwMViMZqamoBwqZXK6jhT2ld0bqdSKdpHbVOgFwlJgV4KLplMdqZeOkbzvTEyVk/T3D25/Hg8zpKly9ncmH5gWV/z9yLlQoFeCi6YronFYqyomNKvSd2GhgaWrt0G7MnfBwN9PieMRaJAgV4GRG/TNd1JJBL8dG36+TUdo3oR6ZoCveRdMFUDfUvX9IYmakW6p+WVkncdqZoOsViMhoaGgnxWZXV8rx8iqVSK91c3FeSzRAYrjeilIIKpmkLKNVHbkb8XkTSN6EVEIk4jesmpvyUBlq7dVrTVL7u2rOucpK2sjvc6X5/db5VEkMFOgV7yoi31GPH4rUA6T86ocd03KJDspZeAJmal7IVK3ZhZvZmtMbMWM9tnPZuZnW9mzZk/L5rZUYFj681shZmlzGx5PjsvpeP91U2dE7CxWIzK6nhR+pFIJDisYS6HNcxlxOiqovRBpNT0OKI3s6HAXcBpQCuwzMwWufvqwGnrgGnu/o6ZnQ4kgWMCx09293/lsd9SgvK5Vl5E8idM6qYOaHH3tQBmNh+YAXQGend/MXD+UmBsPjsp0lf9zdeLREGY1M0YYENguzWzryv/AfwpsO3AE2b2ipl1WZDEzBJmttzMlm/dujVEt0S6V1kd70zf7NqyTuvrpWyFGdFbjn2e80Szk0kH+hMDu09w901mNhr4s5m97u5L9nlD9yTplA+1tbU5319KS1vqsc7guWvLOjj8oOJ2KMvIWH3nCF6lEqSchRnRtwLBJRRjgU3ZJ5lZDXAPMMPdO+9YcfdNmb+3AAtJp4IkAt5f3dS5smXE6KqC3f0qIv0TZkS/DJhoZlXARmAWsNf/0WY2HngY+La7vxHYXwkMcfe2zOvpwE356rwU34jRVRzWMBeARKK015sH8/WgnL2Ujx4DvbvvNrMrgceBocC97r7KzGZnjs8DbgQOAv7LzAB2u3stcCiwMLNvGNDo7o8V5JuIdCN7uafW2Es5CXXDlLsvBhZn7ZsXeH0pcGmOdmuBo7L3S+kp9eWQ/e1fMF8PvcvZ68HhMtip1o2ISMSpBIKElkwm2dy4p1Lkri3rdPepyCCgQC+hNTY27hXcR4yuKlqpg3zQzVRSLhTopVeCq2wGs+APKE3MStQp0EtZ0s1UUk4U6KVHHatONuvJTSKDklbdiIhEnEb0ImhiVqJNgV66lV24LIrLKTUxK1GnQC/d6ihcNmJ01aBfTtkVTcxK1CnQS4/CLqks9TIKIuVKk7EiIhGnQC8iEnFK3YhkUd16iRoFepEA1a2XKFKgLzOqrd69/tStFylVCvSyj3JYOy9STkJNxppZvZmtMbMWM9tniGNm55tZc+bPi2Z2VNi2UnqyH/odxbXzIuWkxxG9mQ0F7gJOA1qBZWa2yN1XB05bB0xz93fM7HQgCRwTsq2UoKiUI86H4ORs8vCNJBKJIvdIpHfCjOjrgBZ3X+vuu4D5wIzgCe7+oru/k9lcCowN21aklFVWxztTV7u2rKOxsbHIPRLpvTA5+jHAhsB2K3BMN+f/B/CnPrYVKSkqjyBRECbQW459nvNEs5NJB/oT+9A2ASQAxo8fH6JbImkDWXph6dptnZ8XdtVSdv+02kkGWpjUTSswLrA9FtiUfZKZ1QD3ADPcfVtv2gK4e9Lda9299pBDDgnTdxERCSFMoF8GTDSzKjMbAcwCFgVPMLPxwMPAt939jd60FRGRwuoxdePuu83sSuBxYChwr7uvMrPZmePzgBuBg4D/MjOA3ZnRec62Bfou0kfJZJLNjXd2bmvtfNe0AkcGo1A3TLn7YmBx1r55gdeXApeGbSulpbGxca/grrXzuWU/oKSxsVGBXgYF3RkrgNbNh9GbFTiqzS+lRGWKRUQiToFeRCTilLoR6aNUKkU8Hu/cbmhoUM5eSpICvUgfVFbHmdK+onM7lUoBKNBLSVKgF+mDkbF6mubuWZIaHNmLlBrl6EVEIk4j+jLVlnqMePxWIJN2GDWu+wYiMmhpRF+m3l/d1JlXjsViukFKJMI0oi9jsViMpqYmQDf45ENwFU5bxRQ9UFxKhgK9SB40NDR0vk6lUrSP2qZALyVDqRuRPEgkEjQ1NdHU1EQsFit2d0T2okAvIhJxCvQiIhGnHL1IAQTr1ldWx5Wvl6JSoC8jwQeM7NqyDg4/qMg9iqaGhgaWrk0/TXPXlnUACvRSVAr0ZST4gJERo6v2Wiki+ZNIJPjp2jFAz3XrRQaCAn2ZCT5gJJH4epF7IyIDIdRkrJnVm9kaM2sxs32GKGY2ycxeMrMPzezarGPrzWyFmaXMbHm+Oi4iIuH0OKI3s6HAXcBpQCuwzMwWufvqwGlvA1cD3+jibU5293/1s68ig1JwYhb0UHEZeGFG9HVAi7uvdfddwHxgRvAEd9/i7suAjwrQR5FBq7I63vnQddjzUHGRgRQmRz8G2BDYbgWO6cVnOPCEmTlwt7snc51kZgkgATB+/PhevL1IachVLyj4QHHQ5KwUR5gRveXY5734jBPc/WjgdOAKMzsp10nunnT3WnevPeSQQ3rx9iIi0p0wgb4VCBYrHwtsCvsB7r4p8/cWYCHpVJCIiAyQMKmbZcBEM6sCNgKzgFALsM2sEhji7m2Z19OBm/raWemdZDK5Vz5YDxgpDcFyxnqguAyEHgO9u+82syuBx4GhwL3uvsrMZmeOzzOzw4DlwCjgEzO7BqgGDgYWmlnHZzW6+2MF+Sayj8bGRlKpVGc1xVgsxoqKKcXtVJkLPlRcDxSXgRLqhil3Xwwszto3L/B6M+mUTrb3gKP600Hpn+DDRUAPGMmnvlzLkbF61pOenG1fq4lZGRiqXikiEnEK9CIiEadAL1JEHROz8XicZDLnLSYi/aaiZiJFoolZGSgK9CJFMjJWT9Pc9PMBOpZbihSCUjciIhGnEX3EBG+SCq6hl9IXvJEKdDOV5I9G9BHTcZMUpNfQ6ylSg0NDQ8NeP5RTqZSqXEreaEQfQdk3SUnpSyQSe43elbOXfNKIXkQk4jSiL2Mqh1B8wf8G6+fu/QxfFT+TfFGgFylBwbkVrbGX/lKgFylBwZy98vXSX8rRi4hEnAJ9BCSTyc56KR2/5ku0qCaO9IdSNxHQ2NjIkqXLGTG6CkaN09r5iOiYqG2rmELHEnvl66UvFOgjYsToKg5rmAtAIvH1Hs6WwUQ1caS/lLoREYm4UIHezOrNbI2ZtZjZPs8/M7NJZvaSmX1oZtf2pq2I9E4wX6+cvYTRY+rGzIYCdwGnAa3AMjNb5O6rA6e9DVwNfKMPbUUkpOz5F+XsJYwwOfo6oMXd1wKY2XxgBtAZrN19C7DFzLKTwz22FZHwVBNH+iJMoB8DbAhstwLHhHz/0G3NLAEkAMaPHx/y7aUnKnMgImECveXY5yHfP3Rbd08CSYDa2tqw71+WgjXnIfPr+6hxxeuQFJVq4khPwkzGtgLBKDIW2BTy/fvTVroQrDkP6bLEldXxovVHiidYx1417KUrYUb0y4CJZlYFbARmAWHvyOlPW+lGds15pWjKk2riSBg9Bnp3321mVwKPA0OBe919lZnNzhyfZ2aHAcuBUcAnZnYNUO3u7+VqW6DvIlL2lMaRXELdGevui4HFWfvmBV5vJp2WCdVWRPJPpY2lKyqBIBIRSuNIVxToRSIqmMYBpXLKmQL9IBFcUplKpTpXWojkojtoJUiBfpDoWFIZi8WIxWIqRSzd0h20EqRAP4hkL6mU8tHdQ8TDCrsiJx+f1ZXsZcD5fn/JTYFepAxoRU55U6AXKQNakVPeFOhFypBW5JQXBXqRMqMVOeVHgb6EaUmlFIJW5JQfBfoSpiWV5aWvhen6skom2Gbz2m3s2rKOivFTAKisjjMyVt+nvvT0WVIcCvQlTksqpdCCJa53bVkHkNdAL8WnQC9S5kbG6jsD++bGOezaso7NjXMASB6+Ubn7CAjz4BERKROV1XFGjK4C0qN7PcgkGjSiF5FO2aN7LcOMBgX6EpLrWbBaaSPZwk5u9ncStLI6zpT2FZ3bWoY5eCnQl5DgKhtAK22kqEbG6mmae2fndjwe1xOsBqlQgd7M6oFfkn4c4D3uPjfruGWOnwF8AFzk7q9mjq0H2oCPgd3uXpu33keQVtlIqVK9nMGrx8lYMxsK3AWcDlQD55lZddZppwMTM38SwH9nHT/Z3WMK8iKDVyKRoKmpiaamJmKxWOfoPh6Pk0wmi9096UaYEX0d0OLuawHMbD4wA1gdOGcG8D/u7sBSMzvAzD7n7m/lvcciUnQa3Q8uYQL9GGBDYLsVOCbEOWOAtwAHnjAzB+5295w/+s0sQfq3AcaPHx+q8+Um7B2QuhNRCi27Gmb26pwVFVN6fdNVd/++C1kjvxyECfSWY5/34pwT3H2TmY0G/mxmr7v7kn1OTv8ASALU1tZmv39kqZ6NDHa5iqS1j9qmu2tLSJhA3wqMC2yPBTaFPcfdO/7eYmYLSaeC9gn05Ur1bGSwy1UkbcnS5Z131+a7do70XphAvwyYaGZVwEZgFpAdjRYBV2by98cA2939LTOrBIa4e1vm9XTgpvx1Pxq00kaipKGhgaVrtwHw4YaVfLhhJe+vbgIU9Iulx0Dv7rvN7ErgcdLLK+9191VmNjtzfB6wmPTSyhbSyysvzjQ/FFiYXn3JMKDR3R/L+7cQkZKRSCT46doxALSlHusM8iqYVjyh1tG7+2LSwTy4b17gtQNX5Gi3Fjiqn30UkUGqu4JpoKJpA0V3xg4wlTmQQsr3iqvs9+vPipdgOWRIp3Uuv/zyzv8f2gIrdUr5ewxGCvQDTGUOpFwFR/eQTut01NJ59tlngWeVyy8QBfoi0OSryN61dJLJJFffnH6dPYELCvz9pUAvIkXX1QQu7Bv4ldfvPQX6AaCbokTCy5Xi6Qjy2Xl9VdAMR4F+APTnpqiuJqVU5kDKRTDwdwT9pWu38eGGlTz77LN7LW7ob+kFKL2J2nz0T4F+gCgvL9J/2UE/+GAUTeh2TYG+QJSuESms7AejdDehW+5BX4G+QFTDRmRgdTWhm2sVT3zprZ2v2/qQ7hlsFOgLSOkakeLIldfPJTvdA9Fc1aNAH0KYyZCe7nhVLXmJmnz8W+3ve4Rpn72KB6Ap8/9gMN0D+67q6SjOBvumf/o7adtd3/M9IaxAnye641Vk8Amme2DfCd4OPaV/Sn2ZpwJ9P+SacFWqRmTwyp7g7Rh195T+yV7m2ZVi/UBQoO8HTbiKlIee0j9hgnz2D4TNgbRQtnzPEyjQ91Jb6jHi8fSvbBrFi0j2E7a6EvYHQvY8QfYPhGDKKCwF+l56f3UTqfc2aBQvZSkqiwW6+x5hv2PY8zomVrPnAw47Nvf5wbt/80WBvgfJZJLNjXtydru2rOPYY2s1iheRgsiVJgpq6mJFTuZJfjkp0OcQ/BUrvc4WPjVuMgAjRldpFC8ig0qoQG9m9cAvST8z9h53n5t13DLHzyD9zNiL3P3VMG2Lpbt8WUdwnzZtGtOmTdunUFIiUVpFj0REutNjoDezocBdwGlAK7DMzBa5++rAaacDEzN/jgH+GzgmZNu86s0MOKSDebZp06bttQwqKnlJESlPYUb0dUBL5kHfmNl8YAYQDNYzgP/JPCR8qZkdYGafAyaEaLuPNWvWEI/He/lV0roL4EHZwVxEJKosHZu7OcHs34F6d780s/1t4Bh3vzJwzqPAXHd/PrP9FPAD0oG+27aB90gAHVH3i8Ca/n21fjsY+FeR+1AqdC320LXYQ9dij1K4Fp9390NyHQgzos81lZv906Grc8K0Te90TwLJEP0ZEGa23N1ri92PUqBrsYeuxR66FnuU+rUIE+hbgXGB7bHAppDnjAjRVkRECmhIiHOWARPNrMrMRgCzgEVZ5ywCLrC0Y4Ht7v5WyLYiIlJAPY7o3X23mV0JPE56ieS97r7KzGZnjs8DFpNeWtlCennlxd21Lcg3yb+SSSOVAF2LPXQt9tC12KOkr0WPk7EiIjK4hUndiIjIIKZALyIScQr0IZjZtWbmZnZwsftSLGZ2q5m9bmbNZrbQzA4odp8GkpnVm9kaM2sxsznF7k+xmNk4M3vGzF4zs1Vm9r1i96nYzGyomf01cz9RSVKg74GZjSNdwuEfxe5Lkf0ZmOzuNcAbwA+L3J8BEyjlcTpQDZxnZtXF7VXR7Aa+7+5fAo4Frijja9Hhe8Brxe5EdxToe/YL4P/SxY1e5cLdn3D33ZnNpaTviSgXnWVA3H0X0FHKo+y4+1sdBQvdvY10gBvTfavoMrOxwNeBe4rdl+4o0HfDzM4GNrr7/xa7LyXmEuBPxe7EABoDbAhst1LGwa2DmU0Avgz8pchdKabbSQ8EPylyP7pV9vXozexJ4LAch/4TuB6YPrA9Kp7uroW7/yFzzn+S/vX9gYHsW5GFLuVRLszsM8BDwDXu/l6x+1MMZnYmsMXdXzGzeJG7062yD/Tufmqu/WY2BagC/jfz5JaxwKtmVufumwewiwOmq2vRwcwuBM4ETvHyugEjTBmQsmFmw0kH+Qfc/eFi96eITgDONrMzgApglJn9zt2/VeR+7UM3TIVkZuuBWncvdoW6osg8QObnwDR331rs/gwkMxtGegL6FGAj6dIeDYPoLu+8yTxk6H7gbXe/psjdKRmZEf217n5mkbuSk3L0EtavgJHAn80sZWbzit2hgZKZhO4o5fEa8P/KMchnnAB8G/hq5t9BKjOilRKmEb2ISMRpRC8iEnEK9CIiEadALyIScQr0IiIRp0AvIhJxCvQiIhGnQC8iEnH/H3GrnN7HIYFGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "from scipy.stats import levy_stable\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "def sample_Zn(alpha, beta, c, n):\n", + " assert n >= 1 and type(n) == int\n", + " \n", + " E_X = alpha*b/(alpha - 1)\n", + " \n", + " samples = [inversion_sample(lambda p: f_inv_pareto(alpha, beta, p)) for _ in range(n)]\n", + " return c*n**(1/3)*(np.mean(samples) - E_X)\n", + "\n", + "alpha = 3/2\n", + "beta = 1\n", + "# TODO: Find proper c in the derivation above.\n", + "c = np.pi/6\n", + "n = 1000\n", + "pdf = lambda x: levy_stable.pdf(x, alpha, beta)\n", + "samples = [inversion_sample(lambda p: sample_Zn(alpha, b, c, n)) for _ in range(1000)]\n", + "compare_plot(samples, pdf, -5, 5, 100)" ] }, {