{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyPPbUR8ohtkP9FJR6sP2FYt"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["# Sampling standard Gaussian random variables"],"metadata":{"id":"856tFswmENXK"}},{"cell_type":"code","execution_count":1,"metadata":{"id":"V0_GOOUdY9X_","executionInfo":{"status":"ok","timestamp":1697483838271,"user_tz":-120,"elapsed":399,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"outputs":[],"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","import math\n","from scipy.stats import norm # To access the CDF of the standard Gaussian distribution"]},{"cell_type":"markdown","source":["Throughout this tutorial, we are only allowed to use Numpy's random number generator to sample uniform random variables."],"metadata":{"id":"OfdXU_IGvjye"}},{"cell_type":"code","source":["def Unif():\n"," return np.random.uniform()"],"metadata":{"id":"65qU4uPjvxW0","executionInfo":{"status":"ok","timestamp":1697483844843,"user_tz":-120,"elapsed":9,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","source":["Unif()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"-LmH8Ve0v0y9","executionInfo":{"status":"ok","timestamp":1697483847416,"user_tz":-120,"elapsed":580,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"90b0b048-a8d1-4877-d291-61733357bbc7"},"execution_count":3,"outputs":[{"output_type":"execute_result","data":{"text/plain":["0.27943576228382394"]},"metadata":{},"execution_count":3}]},{"cell_type":"markdown","source":["The purpose of the tutorial is merely to implement the Box--Muller and Ziggurat algorithms, we will not try to optimise their efficiency."],"metadata":{"id":"x9GuSRDixM9M"}},{"cell_type":"markdown","source":["## 1. The Box--Muller method\n","\n","We recall that the Box--Muller method is based on the remark that if $(U_1, U_2)$ are independent $\\mathcal{U}[0,1]$ variables, then\n","$$X_1 = \\sqrt{-2\\log(U_1)} \\cos(2\\pi U_2), \\qquad X_2 = \\sqrt{-2\\log(U_1)} \\sin(2\\pi U_2),$$\n","are independent $\\mathcal{N}(0,1)$ variables.\n","\n","Complete the function BoxMuller() so that it returns the variable $X_1$ as above."],"metadata":{"id":"e9XL05DXL_CB"}},{"cell_type":"code","source":["def BoxMuller():\n"," U1 = Unif()\n"," U2 = Unif()\n"," return (-2*np.log(U1))**.5*np.sin(2*math.pi*U2)"],"metadata":{"id":"Mklav__ovfcr","executionInfo":{"status":"ok","timestamp":1697483869727,"user_tz":-120,"elapsed":827,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":4,"outputs":[]},{"cell_type":"code","source":["BoxMuller()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"hfv2ACrNwcgN","executionInfo":{"status":"ok","timestamp":1697483872554,"user_tz":-120,"elapsed":259,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"ec79fc0e-20f8-4a66-d7f9-6da70cb86d34"},"execution_count":5,"outputs":[{"output_type":"execute_result","data":{"text/plain":["0.25125458138828327"]},"metadata":{},"execution_count":5}]},{"cell_type":"markdown","source":["We check the validity of the function BoxMuller() by plotting the histogram of 1000 realisations together with the density of the $\\mathcal{N}(0,1)$ law."],"metadata":{"id":"vpE-9zbFxjH8"}},{"cell_type":"code","source":["# Sample\n","N = 5000\n","Sample = np.zeros(N)\n","for k in range(N):\n"," Sample[k] = BoxMuller()\n","\n","plt.hist(Sample, bins=50, density=True)\n","\n","# Theoretical density\n","x_max = max(abs(Sample))\n","xx = np.linspace(-x_max, x_max, 500)\n","den = np.exp(-xx**2/2)/(2*math.pi)**.5\n","\n","plt.plot(xx,den,color='red')\n","\n","plt.show()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":430},"id":"C9FKpNdbx4To","executionInfo":{"status":"ok","timestamp":1697487029771,"user_tz":-120,"elapsed":3163,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"5e8cb37e-e9d7-4af2-bed3-28e51b0d6af5"},"execution_count":40,"outputs":[{"output_type":"display_data","data":{"text/plain":["