{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Random numbers from an arbitrary distribution\n", "K. Reygers, 2020, \n", "inspired by https://tmramalho.github.io/blog/2013/12/16/how-to-do-inverse-transformation-sampling-in-scipy-and-numpy/" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.interpolate as interpolate\n", "import scipy.integrate as integrate" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "# Bose-type distribution\n", "def f(x):\n", " # return 0 if x = 0\n", " num = x**3\n", " den = np.exp(x) - 1\n", " return np.divide(num, den, out=np.zeros_like(num), where=den!=0)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "def get_random(f, xmin, xmax, n_samples):\n", " \"\"\"Generate n_samples random numbers within range [xmin, xmax] \n", " from arbitrary continuous function f\n", " using inverse transform sampling \n", " \"\"\"\n", " \n", " # number of points for which we evaluate F(x)\n", " nbins = 10000\n", " \n", " # indefinite integral F(x), normalize to unity at xmax \n", " x = np.linspace(xmin, xmax, nbins+1)\n", " F = integrate.cumtrapz(f(x), x, initial=0)\n", " F = F / F[-1]\n", " \n", " # interpolate F^{-1} and evaluate it for \n", " # uniformly distributed rv's in [0,1[\n", " inv_F = interpolate.interp1d(F, x, kind=\"quadratic\")\n", " r = np.random.rand(n_samples)\n", " return inv_F(r)" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApv0lEQVR4nO3dd3xUVfrH8c+TSgAhdDWUBAGRoiChRkREmihFkQXERWUXC2DdXQOoKLaw7v6sWLCvqKig0psEgxSl9xpIBAICCqGZPuf3xwwYYpAJZOZMed6vV14zc+feyfcafHJy7rnniDEGpZRSgSvEdgCllFKepYVeKaUCnBZ6pZQKcFrolVIqwGmhV0qpABdmO0BRVatWNbGxsbZjKKWUX1m1atUvxphqxb3nc4U+NjaWlStX2o6hlFJ+RUR+Ott72nWjlFIBTgu9UkoFOC30SikV4LTQK6VUgNNCr5RSAU4LvVJKBTgt9EopFeC00CulVIDTQq+UUgHO5+6MVb4jISmZjMysP2yPiY5iSeL1FhIppc6HFnp1hoSkZH7OPEEj+Yku5Q8xpkdFyM+F0DC46FKoUo/L39hnO6ZSqgS00Aepwq31mOgoljzWEdJSSDw5jpsrboGco5AHLPjjsRsiQ+F/H0LjPtD4FihTwavZlVIlo4U+SGVkZpGe1AOAwaOehzcfh4ObSAgpD436QN3r4OKmULEWhJWBghw4vh8ObuH9iZ9yb+ZmmP4gzHsS4u+ENsPgohpWz0kpVTwt9MHsxEGY/iAfRcyCvDjo/SZtJ0WxrVefP+4bEgWV60LluiTlw70jboSMVbBsPCx9DVa8D9c+Cm3uh7BI75+LUuqsdNRNkGofsh7eaAupC3gubyAMWw7NBpJDxDmPjYmOInbkLGJfP0DCzjtg+EqIaw/fPgVvXQMZqz1/Akopt2mhD0bL3+HD8HFQvgbck8I7BTdB2LkL/ClLEq8nPakH6Uk9nP38VS6DAZ/B7VMg5wS8ewMsehEcDg+ehFLKXdp1E0yMgYXPw6J/s9DRnBuGzIDI8sCu8/7ImOgoYhNnnn5dgad4JvxDeiU/CxlroM9berFWKcu00AeLQkWe5oMYuqwbuyLLA2cW65joqBJ9bLHj6U0/xjz+IE9vnwjvduK2Yw+z4ljF05+vY/CV8i4t9AHu1DDKu0JnMyb8Y2g+CG5+Dcey2af3KfXCK8JHBV15+q+3wRd3MD5nJNUfmgEXNz2j9a+U8g7tow9wGZlZpN8pjAmfyJyClnDzaxDi+R97THQUsW8fo1PmKIyEwgc3wk9LPf59lVJ/pC36ANdQdsOUsXBpcx7adT9bvVDkochfCUdvhI/7wMS+xMujQA+vZFBKOWmLPpDlnGB8+CsQWQEGTCIbS+PbK9aEwTOgwiV8EPEi7F1pJ4dSQUoLfaAyBmY+Qqz8DH3fs3/X6kU1YPB0fjUV4ONbYN9au3mUCiJa6APV+s9h/ee8mn8LxF4D/D66JjZxZolH15SKCpcyMHe0c7jlJ7fBkZ+8n0GpIKR99IHo+AGY/RjUas1rO/rwsGuzLwxr3EdVGDQF3usMn/SFu+dC2cq2YykV0LRFH4hm/QPysqDn6zh88Udc7XLo/xkcSYdJAyEv23YipQKaD1YBdUE2T4Mt0+C6x6BaA9tpzi42AXq/CbuXwbQRzmsKSimP0K6bQJL7G8wZCTWaQrsHbKcp1pl34VZmScfRsPA5uLQZtB1mN5xSAcqtFr2IdBORbSKSKiKJxbz/iIhsFpH1IrJAROoUem+wiOxwfQ0uzfCqiKWvwrG90H0chIbbTlOsP0yI1v4f0PAmmPcE7PrOdjylAtI5W/QiEgqMBzoDe4EVIjLNGLO50G5rgHhjzG8ich/wb+AvIlIZGAPEAwZY5Tr2SGmfSNA7uhcWvwyNepMwKYeMzPObu8abYqKjiB01m3L0YUbUOuK+vAuGfgeV6pzzWKWU+9zpumkFpBpjdgGIyCSgF3C60BtjFhba/wdgkOt5V2C+Meaw69j5QDfgswuPrs4wfwxgoPNYMsZtPL16lC8rPAqow8ijpEQ+DZ/fDkO+hfAyFpMpFVjc6bqJAfYUer3Xte1shgCnZsxy61gRGSoiK0Vk5aFDh9yIpM6wbw1snMxr2d2IHbfRp1vxZ/OTuRhueQd+3gDzHrcdR6mAUqoXY0VkEM5umg4lOc4YMwGYABAfH6/DL0pq4fNkmnKMePxVRpSpaDvN+WvQFdoOh2WvO1esatTLdiKlAoI7LfoMoFah1zVd284gIjcAo4Gexpickhyrzk9CUjK3jHwJdsxjUngf8Ocif0qnMRDTAqaOcI6zV0pdMHcK/QqgvojEiUgE0B+YVngHEWkOvI2zyB8s9NZcoIuIVBKRSkAX1zZVCjIys/iq4UIoW5V7/znOdpzSERYBfd93Pp98N+Tn2s2jVAA4Z6E3xuQDw3EW6C3AF8aYTSIyVkR6unZ7ESgPfCkia0VkmuvYw8AzOH9ZrADGnrowqy5ca9kCaSlwzcOuJQEDRKVY6PkqZKyChc/aTqOU33Orj94YMwuYVWTbk4We3/Anx74PvH++AdXZDQv7BspVh5ZDbEcpfY17w87BsORVaNAN6rSznUgpv6VTIPir/eu5NnQDtLkXwv1vlM2fSUhKJjZxJo2WduAnRzX2fTgYso/ZjqWU39IpEPzVklc4YcpQPv5u20lKRdEFyk/fB7C7JgXvdYW5I6HXeIsJlfJfWuj90ZF02PQ1nxZ0Y2hUJdtpSsVZp1Cu3Zo3C3oyfM1EuPxGaOj7N4Ip5Wu068YfLRsPEsL7+d1sJ/GKV/JvhYuvhGkPwImD5z5AKXUGLfR+pusL0/ntx4+YnNeW0OiatuN4RR5hzrtmc47DzEdtx1HK72ih9zPtTsylrOTQ975nfGLFKG+IiY4i9v92kpTdxznX/qZvbEdSyq9oH70/cTi4I3Q+1GzpnL89SJz+hVbQlfVP/8iVs/4BcdfqEoRKuUlb9P5k10LqhvwMrYbaTmJHaBj/yrsHso7AnD8si6CUOgst9P5kxbv8YioE9WRfW01taP8orP8ctutsGkq5Qwu9v8jcDdvnMKmgI4RF2k5jV/t/QPVGMP0hyD5qO41SPk8Lvb9Y+QEAn+SfdbaJoBATHUXs4/PpuWcABcd/di5BqJT6U3ox1h8U5MPaT6B+V/avr2I7jVWFRxq9/fgP3LP6I2ja13lx1iUhKdm5Hi3OXwzBMjpJqbPRFr0/SP0WThyAq++wncSnTCp3B+mOGuz8YCgdX5hzentGZtaZC5ArFeS00PuDNR9DuWpQv4vtJD5l4cjuxP71TS4L2U+vk5/bjqOUz9JC7+tOHILtc+DKv0BouO00vqdeJ2jSl/tCp8EvO2ynUconaaH3dRu+AEc+NB9kO4nv6vo82UTAjIfB6JLDShWlhd6XGQOrP4aYeKh+he00vuuiGiTlD4D072HdZ7bTKOVztND7sn2r4dAWaH677SQ+b1JBR6jVGuaOphK6SIlShWmh92Xrv4DQSGh8i+0kPs8QAje9BDnHGBmmrXqlCtNx9L6qIB82fgUNupLwyuozxoWrs6jRGNqNoN/ilyDte4hrbzuRUj5BC72PeijpVV7OO8i96+pCBX5fWk8V69RShGW4kgVRNYiZ8TDct+QPSxTqzVMqGGmh91Htsr6DchV4a3RiwC3+7QlnFPAdF8EnfWHZ6yxJ/H2hklMFX6lgo330vigvm26hy+GKm7XIn4/6naHhTZDyonMyOKWCnBZ6X7RjHhUkyzmHizo/3ZKcj3NG2s2hlA/QQu+LNnzJIVMRYq89976qeNG1oMM/YesM2DHfdhqlrNJC72uyj8H2ucwoaAOhegnlgrQdDlXqwax/Ql627TRKWaOF3tdsnQEFOUwraGc7if8Li4QbX4QjabD0VdtplLJGC72v2TwNKtZijalnO0lguOx6aNQbvv8vNeWg7TRKWaGF3pfkHIedyc7RNojtNIGj6/MgoYwJ+5/tJEpZoYXel+yYBwU5rkKvSk3FGLjuMTqHroZts22nUcrrtND7ks3ToFx15+RcqnS1vo8djhiY/Rjk6apTKrhoofcVeVnOYYANe0BIqO00gScsgifz74TMn2DxS7bTKOVVWuh9xc5kyDsJjXraThKwljkaQ5O+sPhl+HWn7ThKeY0Wel+xZTqUiYZYnXHRU2Kio2i1siMn8oVl4/+mq1GpoKF35PiC/FzYNotZuVdz/+h5gE5H7AmnJz5bmknbeaNh60y44ia7oZTyAi30viD9e8g+yle5V+t0xN7Q+h62znmLhnMSnePsI8raTqSUR2nXjS/YMg3Cy/G9o6ntJMEhNJwn8u6Co3vg+//YTqOUx7lV6EWkm4hsE5FUEUks5v1rRWS1iOSLSN8i7xWIyFrX17TSCh4wHA7n2O76nckhwnaaoLHCNIQr+8OSV+GXVNtxlPKocxZ6EQkFxgPdgUbAABFpVGS33cCdwKfFfESWMaaZ60uHlBS1fw2cOACXd7edJKjEREfRcvm1HCsIY/kbd+uFWRXQ3OmjbwWkGmN2AYjIJKAXsPnUDsaYdNd7Dg9kDGzb54KEQL3OwA+20wSN0xdmfzxBq9n/gs1ToXFvq5mU8hR3um5igD2FXu91bXNXGRFZKSI/iEjv4nYQkaGufVYeOnSoBB/t/7Yt+pIVBfWIfeYHHWljQ/wQNjnqwNxRkHPCdhqlPMIbF2PrGGPigYHAyyJyWdEdjDETjDHxxpj4atWqeSGSjziaweVmFy27DCQ9qYcuXG1DaJjzwuyxDFj0b9tplPIIdwp9BlCr0Ouarm1uMcZkuB53Ad8BzUuQL7DtmOt8bNDNbo4gt9o0gGaDYNl4OLTNdhylSp07hX4FUF9E4kQkAugPuDV6RkQqiUik63lVIIFCfftBb9scdjuqQbWGtpMEtZjoKK7+4RqOFkSw9NU7iU2cQWziTBKSkm1HU6pUnPNirDEmX0SGA3OBUOB9Y8wmERkLrDTGTBORlsDXQCXgZhF52hjTGLgCeNt1kTYESDLGaKEHyP0N0lJY4OjAXaJzz9t0ustsxUnazXyU9H7Z0LQvsYkz7QZTqpS4dWesMWYWMKvIticLPV+Bs0un6HFLAb0LqIiEpGQaHlvMexHZrItqYzuOOqXFXbD6Y5g7Gup3sZ1GqVKjd8ZakJGZxXttf4WI8rz82HDbcdQpIaHQ4/+c9zWkjLOdRqlSo3PdWGGc4+cvux7C9G5Yn1KzBbQYDD+8yeXyhz9SlfJLWugtaCzpcHy/3g3rqzqNgc3TGMf/iE2sBQgx0VE6/FX5Le26seD6kDWAuO6GVT6nbGW44Smamc2kDzhBelIPMjJ1+UHlv7TQW9AhdD1c2hzKB9HNYf6m+R0Q0wLmPQ5ZmbbTKHVBtNB7W9YRmssOqNfJdhL1Z0JCoMd/4eQv8N0LttModUG0j97bdqUQKgbq3WA7iTqXS5tDyyGwfAKNpLbtNEqdN23Re1vqtxwzZSEm3nYS5Y7rH4eoyowN/9C5doBSfkgLvTcZA6kLWOxoAqH6x5RfiKoEnccSH7Id1hW33IJSvk8LvTcd2grH95HiuMp2ElUSVw1gpaMBzH8Sso7YTqNUiWmh95KEpGSefeU1ALaXb2U5jSqRkBDnVMZZR2DBM7bTKFViWui9JCMzi8cvz4BqDfl6ZD/bcVQJbTF1oNVQWPk+7FtjO45SJaKF3kuiyIaflupoG3/WcRSUrw4zHgZHge00SrlNC72XtA7ZAgW5On7en5WpCF2fd7bol79jO41SbtNC7yUdQtZDWBTUbmc7ijoPMdFRxCbOJPaTMiwLaQ7Jz8BRtxdaU8oqHePnJdeGrIfYayC8jO0o6jwUntDsmpGHWFxuJMz+FwnpQ07Pg6MTnylfpS16bzicxmUh+7V/PkCYinV4Ias3bJ1B+4IfSU/qoROfKZ+mLXpv2Olae1T75wPCksTroaA9TNhA0m//g5zhEHmR7VhKnZW26L0hLYV9pjJUqWc7iSotoeFw8yvOdQWSn7WdRqk/pYXe0xwOSPuepY4moIuAB5aa8dDyb/Dj25CxynYapc5KC72nHdgIWYdZWtDIdhLlCZ2egPI1YPqDhKJj65Vv0kLvaWkpACxxNLEcRHlEmYpw47/h5w3cFTrHdhqliqWF3tPSFkGV+hygsu0kylOu6AkNuvFI2GTI3G07jVJ/oIXeg659YR4nti/i4wN1iImOsh1HeYoI3PgfEGHhfwcRmziDhKRk26mUOk0LvQdVPbaJ8pLNHQP/qjfSBLroWpTt9jQdQ9eRPuCEjqlXPkULvQclhGwEBGLb246ivKHVUKjVGmY/RjUybadR6jQt9B6UELoJLm4KZbV/PiiEhEDP1yEvi6fDP7SdRqnTtNB7Su5vNJcdULeD7STKm6o1gOsSuTF0OWz6xnYapQAt9J6z5wciJR/itNAHnXYPsN4RB7P+Ab8dtp1GKS30HrMrhTwTCrXb2k6ivC00jMfyhjqXHpyTaDuNUlroPSZtEWtMPYgsbzuJsmCLqQPtH4X1n8P2ubbjqCCnhd4TsjJh/1qWORrbTqIsiYmOov7cxmxz1OTgp/dB9lHbkVQQ00LvCT8tAeNgSYEW+mC1JPF6diT14vKh/6OKOQzznrAdSQUxLfSesCsFwqJYY+rbTqJsq9mCdwpugtUfwY75ttOoIKWF3hPSUqBOW/J0XRcFvJR/K1RvBFOH6ygcZYUW+tJ2/AAc2qrDKtVpOURAn7fgt19g5qOntyckJTsXHE+cqXPjKI/SQl/a0hY5H+OutZtD+ZZLroLrEmHTV7BhMgAZmVm63qzyCrcKvYh0E5FtIpIqIn8YGCwi14rIahHJF5G+Rd4bLCI7XF+DSyu4z0pLcc5RfslVtpMoX5PwMJukAZmTH6BV4kSd0VR5zTkLvYiEAuOB7kAjYICIFF0uaTdwJ/BpkWMrA2OA1kArYIyIVLrw2D4sLcU5iVlIqO0kyteEhjE8+x6iwx0sb/INSx7raDuRChLutOhbAanGmF3GmFxgEtCr8A7GmHRjzHrAUeTYrsB8Y8xhY8wRYD7QrRRy+6bDaZC5myc3VCE2caa22BTgHFN/qi8+t2Jd6DwWUufDqg9tR1NBwp1hITHAnkKv9+JsobujuGNjiu4kIkOBoQC1a9d286N9kKt/fuxDwxhb7XLLYZSv+MNaBI7rYNtMmDvaOeld5bpWcqng4RMXY40xE4wx8caY+GrVqtmOc/7SUjhgoqFqA9tJlC8LCYFe4yEkDL66BwrybSdSAc6dQp8B1Cr0uqZrmzsu5Fj/YgykLWKpo7FzaTml/kzFmnDzS7B3OaQk2U6jApw7hX4FUF9E4kQkAugPTHPz8+cCXUSkkusibBfXtsBzcAucPOQs9Eq5o8mt0GwQLPoPbUI2206jAtg5C70xJh8YjrNAbwG+MMZsEpGxItITQERaishe4DbgbRHZ5Dr2MPAMzl8WK4Cxrm0BJSEpmadffQOAneVbWE6j/Er3cVDlMl4Kf0PvmlUeI8YY2xnOEB8fb1auXGk7RonEJs4k/aqJzlb9g2ttx1H+Zt9act++nogrusNfJmrXnzovIrLKGBNf3Hs+cTHW34VSAOmLddlAdX4ubca/8/vD1hmw8n3baVQA0kJfCppIGuQc02kP1Hl7r6A7XHY9zB3l/MtQqVKkhb4UJIRscj7RiczUeTKEQO+3IPIimDwE8nTuG1V6tNCXgnYhG6FGEyhX1XYU5adioqOIfW4lg4/cDQc3ORcWV6qU6ITpFyovm/iQ7RD3d9tJlB/7/e7ZHrz6+DYeWDMRareD5rdbzaUCg7boL9Te5ZSRPO22UaXmy3KDWFLQmKxvHuKvz+vFWXXhtNBfqF0p5JsQqNPOdhIVIL4f2ZmEx74h6qLKPJU9DrKP2Y6k/JwW+guVtoj1pi6UqWA7iQok5atD3/epLQdh2gjnFBtKnSct9Bci+xhkrGKJo4ntJCoQxSbwYv5fYPM38OPbttMoP6aF/kLsXgamQOe3UR4zs3xf5he0IHf2KO55/g3bcZSf0kJ/IXalQFgZVjvq206iAtTikTfQedQUIqrU4ZmccXBsn+1Iyg9pob8QaSlQqzU5RNhOogJZVCXo/yllyYZJt0Netu1Eys9ooT9fJ3+BAxt12gPlHdWv4NmIh2DfaiaPvY2EFxbYTqT8iBb68/TEy28C0Ht2uK4Nq7wiafRI6JBI39BFdDnxte04yo/onbHnqWHWGihXgW+eGAah+p9ReUmHx+DARkZv/YSBo2qx1NGEmOioP65Lq1Qh2qI/T+1CNkKdBC3yyrtCQqDPW4RVa8CnFd8k/Z8NycjUCdDUn9NCfz4y9xAXckDnn1d2RF4E/T8FBD7tR0VO2E6kfJwW+vORtsj5qBdilS1VLnMW+8zdvB3xEuTn2k6kfJgW+vORlsIhUwGqN7KdRAWzOm2h1xu0CdkC0x/QaRLUWWmhLyljIG0RPzga6dqeyr4rb+P/8vrCus9g0X9sp1E+Sgt9Sf2yHY7vZ7Gjqe0kSgEwpfxAphRcAwufZcxzT9uOo3yQDhkpqV0pACzR+W2Uj1gyshPkXwMf92FU+quw6zodKKDOoC36ktr1HVSKZa+pbjuJUr8Li4T+n5BuLnZOk7Bvre1EyodooS+JgnxIX6yrSSnfFFWJf0U9xd6cSA69fTP9nv/EdiLlI7TQl8T+tZBzVP8sVj5r6sjbqDliDtXKhfOf7DFw/GfbkZQP0D76ktj1nfMxrgPwo80kSp1d1fpw+5dUmdAdJt4Kd84k4ZXVp++g1SkTgo+26Etg1Xdfs9lRh9hnftSJzJRvi2nBPXmPwKFt8Fl/DmceIT2pB+lJPXTKhCCkLXp35WXRpGAbkQn3kN61h+00Sp3TYkdTuGUCTBnCu+HHIO9GCNcGSjDSFr27dv9ApORB3etsJ1HKfU1ugd5v0jZkM3w+CPJzbCdSFmihd9eu78gzoVC7re0kSrklJjqK2MSZxH52EePC74fUb+GLwYSTbzua8jLtunFXWgqrTX1aR5a3nUQpt5x5wbUHrKgLMx/lnbJHqJcI+YTphdkgoS16d/x2GPatZWmB3g2r/FjLv0G3JK5z/EDq1V+R/lwXvTAbJLRF7470xYBhiaMxD9vOotSFaHMfOPJh3uNQkEsk/WwnUl6ghd4du76DiPKsza5nO4lSF67dCAgrA7P+wTvheyC3K0SUtZ1KeZB23bgjLQXqtCNffy+qQNHq79BrPAkhG+GTvpBz3HYi5UFauc7l6F74NRXi74YNtsMoVYqaD+LpWTt5Mv1lNj7Xgb/mPsYxnIMN9CJtYNEW/bm4piXWicxUIBo7egxhAybSLHwP62u/TPrIZnr3bAByq9CLSDcR2SYiqSKSWMz7kSLyuev9H0Uk1rU9VkSyRGSt6+utUs7veTuToVx1XTZQBa6GPeD2LyFzD7zXGQ5usZ1IlbJzFnoRCQXGA92BRsAAESla9YYAR4wx9YCXgHGF3ttpjGnm+rq3lHJ7h6OAzI1zmXKsAbGjZuv8Nipw1b0O7prlHJHzfldaylbbiVQpcqdF3wpINcbsMsbkApOAXkX26QV85Ho+GegkEgALqu5fSzTHufW2O0lP6qF9liqwXXIlDJkP5aozMeIF2PSN7USqlLhT6GOAPYVe73VtK3YfY0w+cBSo4novTkTWiEiKiLS/wLzelboAhxG4rKPtJEp5R6U6MGQeG00sfHknLH0djLGdSl0gT4+62Q/UNsb8KiItgG9EpLEx5ljhnURkKDAUoHbt2h6OVAKpC9hoYrmyXFXbSZTynrKVuT13FFubfwXzRjN9QTKPnBxMnk6Z4LfcadFnALUKva7p2lbsPiISBlQEfjXG5BhjfgUwxqwCdgINin4DY8wEY0y8MSa+WrVqJT8LT8jKhL0rSHFcZTuJUl5XJTqauDV/4ZX8PtxcsIAdDd4i/YnWOhrHT7nTol8B1BeROJwFvT8wsMg+04DBwDKgL5BsjDEiUg04bIwpEJG6QH1gV6ml96S0FDAFLCq4khG2syjlZb+32m+GDTfB1GEwoSOXy/1Wc6nzc84WvavPfTgwF9gCfGGM2SQiY0Wkp2u394AqIpIKPAKcGoJ5LbBeRNbivEh7rzHmcCmfg2ekLoDICqwxOu2BCnJN+zpH5BTkMiXiKdg81XYiVUJifOxCS3x8vFm5cqXdEMbAy03hkquIXTuQ9CRdUUopju1nzX9uonlIKrQZBp2fhtBw26mUi4isMsbEF/ee3hlbnF+2w9E9UK+T7SRK+Y4Kl9Av90lodQ/8MB4+7AHH9tlOpdyghb44qQucj5dpoVeqsOrRFxG7qAMjcofz25518FZ72LnQdix1DjqpWTGWzf+C6o5L6DRuo94Nq1Qhv1+k7UGnkXVYUPVd+LgPtH8ErhupXTk+Sgt9UbknubpgI5Ht/k56N+2bV+pssivWo9Hex3gq7CP6ff9f1qV8zUN5w8itWFfH2vsYLfRF7VxIpORBg262kyjl034v5rfC5qlcNe0BFhY8QeLxgWA6QgDMghIotNAXtX02x0xZKtRpZzuJUv6jUS+IiYdv7iUp7V3mPbmG0XlDOET06V30rlp79GJsYQ4HbJ9HiuNK7WtUqqQqxsAdU6HLs3SJ2MCK6FGkDzhO+gs36hz3lmmhL2zfGjh5kG8LrradRCn/FBLiXJP23sVQtQF8fQ982s+5UpuyRgt9Ydtng4TwnaOZ7SRK+bdql8Pdc6DrC5C+GMa34fbQb51/NSuv00Jf2LY5UKsNR13rZiqlLkBIKLS9H+5bCjHNeS78fdY+Fc/NI18jISnZdrqgooX+lMw9cGADXK6jbZQqVZXj4K/ToM8Eml10nOmRT3DfidfhN/+Y9ioQaKE/Zfsc52OD7nZzKBWIROCqv8CIldD6XgaEJsNrLWDVh+AosJ0u4OnwSpcVcz6miuNirv/vdmKiy9qOo1RgKlMRuidx95r6DDv5Nq2mP8i2qS/yQv4A17Ux0WGYHqCFHuDkrzQv2EBYh4dJ73ST7TRKBbyPRg0Bczdsnsrl3z7Fh0dehLhrofMzxL5adF0jdaG06wZg6wzCxOG86UMp5R0i0Lg3DFsO3cbBzxthQgdeCh8Pv+60nS6gaKEH2DyVdEcNuPhK20mUCj5hEdDmXnhgDSQ8RLeQFfB6S/j6Pi34pUQL/W+HIS2F2Y5WOjeHUjZFRUPnp+lX5i3eyetK9tovyX8tHr65Hw77xwqkvkr76LfNAkc+swpac5/tLEoppo+8FbgVjh/g3XEP8reNU2DdJOZLW17O6s4mE6cXbEtIW/Sbp0J0bTaYONtJlFKFXVSDZ/PvgAfXQdv7aVOwipmRo0lv+Cb1jy1zLvmp3BLchf63w87VcRr1ArTbRilfExMdRexzq4hNbsutZd6Bzs/Arzv5MOLf8EZbWPUR5J60HdPnBXfXzcYp4MiDpv0geY/tNEqpIortnml9Lw89OYaXQ76H6Q/AvCfgqv7Qcohzjh2XhKTk0zNmBntXT1C36DfNnsAWRy1iX9mtSwYq5S/CIvjGcQ3c+z3cNQcadIFVH8D4VvBBD2cDLi+bjMws0pN66BTJBHOL/pcdNDbboeszpCfoTVJK+ZOY6ChiR85yvepDFa7nttAU7ti9kJif7oYyFXkuLB72VIOaLa1m9QXBW+jXf06BEUKb3mY7iVKqhIrvhhlIXOJ00v5WFtZN4pb138B7C6ByXUaEXg1HmkClOt6O6hOCs+vG4YB1n7PY0RQqXGI7jVKqlFwaXY7Yd7OJXdGbnmU+hF7joUIMj4ZPhleuhAkdYfHLcDjNdlSvCs4W/a6FcHQ3UwpupoPtLEqpUlNsS7/5IG59fhLxJ7+j+94fabZvDHw7xnknfOPecEUvqFrP61m9KTgL/Yr3oGxV5mS3sp1EKeUFU0b1B/oDcOvzk2h2chE37ltOi5/HwoKxpDlqsCK8Jf0G3A11EiC8jN3ApSz4Cn3mHueSgQkPkfutLgCuVLApXPQ5mgHbZhG3Yx6XbJ8LE2dAeFmI6wD1O8NlHaFSnN9PjxJ8hX7VB87H+Lvg2w12syil7KoYA63+Dq3+zlWJX7Pt7nKwYx77Vk7l0u2zAfiZqlx8VWeIbQ9x7SG6tuXQJRdchT4vy3knXYNurh+WFnqllFMOEc4x+Q260O7760h/tAGkpbBq+hf02DEP1n3m3DG6DsS154UNFUk+WYdUcymXRpfz6RuygqvQr/4YfvuFfutbsHzdTL1JSil1Wkx0FLGJM13Py0K1BlCtAc8nX8bwwydpIHvpVm4HD198ALbOZGT+EUZGApEVWHSyDiQvg1qtIKYFlK1s92SKEONjEwPFx8eblStXlv4H5+fCq81ZkVmWlk8v9/s+N6WU9/0+rYKhTYUjTLoxFPauYMeqZOo6fiJUXPW0Ym245Eq4uKlzdM8lV0KFGI/WHRFZZYyJL+694GnRr/8cju1lfP4/+VCLvFLqPBTbPdNsIPVvAnJOwL7VvPnJ59Q8nEqjI6uI2zqTEFzFP6qys/BXb8S4VbDiZDVSTQzloqt7vNsnOFr0uSedK85fdAmxux4lPUmnPFBKed4NL8ykwrHtNA5Jp7Gk0yjkJ+pLBlGSe3qfQ6YC1eKugqoNnN0+zW8/r+8V9C36d//9CH/L30/fX/7u7HtTSikv+HZkD6DHmRsdDji6Bw5tg1+28cOChVy6azf101axc9Uymp9nof8zgV/oD6cxMO9raNSTyX951HYapVSwCwlxzrlTqQ406MLN7UY4txvDoJFT2OSJb+mBz/QdDgdMHUY+IdDtBdtplFLq7EQ4iWdGAgZ2oU8eCz8t4Zn8O6BiTdtplFLKCrcKvYh0E5FtIpIqIonFvB8pIp+73v9RRGILvTfStX2biHQtxex/7se3YfFL0OIuvizQqcuUUsHrnIVeREKB8UB3oBEwQEQaFdltCHDEGFMPeAkY5zq2Ec5JJRoD3YA3XJ/nGcbAkXSYOgxm/4t5BS2ot6SjXoBVSgU1dy7GtgJSjTG7AERkEtAL2Fxon17AU67nk4HXRURc2ycZY3KANBFJdX3estKJX8ixfTC+DeQcBQnljfye3P/0B6SGBv71ZqWU+jPuVMEYoPDK2XuB1mfbxxiTLyJHgSqu7T8UOTam6DcQkaHAUNfLEyKyza30xasK/AITGfbcxAv4GL/iOuegEWznC3rOwaKqjDvvcz7r8lk+0dw1xkwAJpTGZ4nIyrPdNBCogu2cg+18Qc85WHjqnN25GJsB1Cr0uqZrW7H7iEgYUBH41c1jlVJKeZA7hX4FUF9E4kQkAufF1WlF9pkGDHY97wskG+fcCtOA/q5ROXFAfWB56URXSinljnN23bj63IcDc4FQ4H1jzCYRGQusNMZMA94DPnZdbD2Ma/kW135f4Lxwmw8MM8YUeOhcTimVLiA/E2znHGznC3rOwcIj5+xzk5oppZQqXYF9Z6xSSikt9EopFegCptCfa5qGQCMitURkoYhsFpFNIvKg7UzeIiKhIrJGRGbYzuINIhItIpNFZKuIbBGRtrYzeZqIPOz6d71RRD4TkTK2M5U2EXlfRA6KyMZC2yqLyHwR2eF6rFQa3ysgCr2b0zQEmnzgUWNMI6ANMCwIzvmUB4EttkN40SvAHGNMQ+AqAvzcRSQGeACIN8Y0wTkIpL/dVB7xIc6pYQpLBBYYY+oDC1yvL1hAFHoKTdNgjMkFTk3TELCMMfuNMatdz4/j/J//D3cdBxoRqYlzJYd3bWfxBhGpCFyLc2QbxphcY0ym1VDeEQZEue7LKQvss5yn1BljFuEcpVhYL+Aj1/OPgN6l8b0CpdAXN01DwBe9U1yzhTYHfrQcxRteBv4FOCzn8JY44BDwgau76l0RKWc7lCcZYzKA/wC7gf3AUWPMPLupvKaGMWa/6/nPQI3S+NBAKfRBS0TKA1OAh4wxx2zn8SQRuQk4aIxZZTuLF4UBVwNvGmOaAycppT/nfZWrX7oXzl9ylwLlRGSQ3VTe57rptFTGvwdKoQ/KqRZEJBxnkf/EGPOV7TxekAD0FJF0nN1z14tIoM9ctxfYa4w59dfaZJyFP5DdAKQZYw4ZY/KAr4B2ljN5ywERuQTA9XiwND40UAq9O9M0BBTXNNDvAVuMMf9nO483GGNGGmNqGmNicf6Mk40xAd3SM8b8DOwRkctdmzpx5hThgWg30EZEyrr+nXciwC9AF1J4OpnBwNTS+FCfmL3yQp1tmgbLsTwtAbgD2CAia13bRhljZtmLpDxkBPCJqxGzC7jLch6PMsb8KCKTgdU4R5etIQCnQxCRz4DrgKoishcYAyQBX4jIEOAnoF+pfC+dAkEppQJboHTdKKWUOgst9EopFeC00CulVIDTQq+UUgFOC71SSgU4LfRKKRXgtNArpVSA+39CwL8lQxbE7QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from scipy import integrate\n", "ran = get_random(f, 0., 10., 100000)\n", "plt.hist(ran, bins = 100, histtype='step', density=True);\n", "x = np.linspace(0, 10, 1000)\n", "norm, norm_err = integrate.quad(f, 0., 10.)\n", "plt.plot(x, 1/norm * f(x));" ] } ], "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 }