{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyMfWViYOFECvs7d148+pMyT"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["# Variance reduction"],"metadata":{"id":"856tFswmENXK"}},{"cell_type":"code","execution_count":66,"metadata":{"id":"V0_GOOUdY9X_","executionInfo":{"status":"ok","timestamp":1698055014098,"user_tz":-120,"elapsed":627,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"outputs":[],"source":["import numpy as np\n","import matplotlib.pyplot as plt"]},{"cell_type":"markdown","source":["## 1. Control variate\n","\n","For $t>0$ and $x \\in R$, we set $$f_t(x) = \\frac{1}{1+t x^2}, \\qquad g_t(x) = 1-t x^2.$$"],"metadata":{"id":"e9XL05DXL_CB"}},{"cell_type":"code","source":["def f(t,x):\n"," return 1/(1+t*x**2)\n","\n","def g(t,x):\n"," return 1-t*x**2"],"metadata":{"id":"Mklav__ovfcr","executionInfo":{"status":"ok","timestamp":1698047197420,"user_tz":-120,"elapsed":233,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":8,"outputs":[]},{"cell_type":"markdown","source":["For $X \\sim \\mathcal{N}(0,1)$, our aim is to estimate the value of $\\mathcal{I}:=E[f_t(X)]$, using $Y=g_t(X)$ as a control variate. Note that $E[Y] = 1-t$."],"metadata":{"id":"aq08HibqKlOW"}},{"cell_type":"markdown","source":["With the notation of the lecture notes, complete the following code so that the function MCandCV(t,X) takes in input a number $t>0$ and an array $X$ which contains $n$ realisations of independent standard Gaussian variables, and returns the list $(\\hat{\\mathcal{I}}_n, \\hat{\\sigma}_n, \\hat{\\mathcal{I}}^\\mathsf{CV}_n, \\hat{\\sigma}^\\mathsf{CV}_n)$."],"metadata":{"id":"0SEbykCDKWh3"}},{"cell_type":"code","source":["def MCandCV(t,X):\n"," n = len(X)\n"," fX = f(t,X)\n"," Y = g(t,X)\n","\n"," mean_fX = np.mean(fX)\n"," mean_Y = np.mean(Y)\n","\n"," var_fX = np.mean((fX-mean_fX)**2)\n"," var_Y = np.mean((Y-mean_Y)**2)\n"," cov_fX_Y = np.mean((fX-mean_fX)*(Y-mean_Y))\n","\n"," beta = cov_fX_Y/var_Y\n","\n"," estim_CV = np.mean(fX-beta*Y)+beta*(1-t)\n"," estim_var_CV = var_fX - cov_fX_Y**2/var_Y\n","\n"," return [mean_fX,var_fX**.5,estim_CV,estim_var_CV**.5]"],"metadata":{"id":"hfv2ACrNwcgN","executionInfo":{"status":"ok","timestamp":1698048832693,"user_tz":-120,"elapsed":246,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}}},"execution_count":19,"outputs":[]},{"cell_type":"markdown","source":["We use this function to plot the estimated values of $\\mathcal{I}$ as a function of $t$ obtained with both the naive MC estimator and the CV estimator, together with respective $95\\%$ confidence intervals."],"metadata":{"id":"RwL3Fk3gnzgU"}},{"cell_type":"code","source":["# Sample\n","nX = 1000\n","X = np.random.normal(0,1,nX)\n","\n","# Values of t\n","nt = 100\n","t = np.linspace(.1,10,nt)\n","\n","# Estimators and standard deviations\n","estim_MC = np.zeros(nt)\n","sigma_MC = np.zeros(nt)\n","estim_CV = np.zeros(nt)\n","sigma_CV = np.zeros(nt)\n","\n","for i in range(nt):\n"," [estim_MC[i],sigma_MC[i],estim_CV[i],sigma_CV[i]] = MCandCV(t[i],X)\n","\n","# Plot of estimators and 95% confidence intervals\n","plt.plot(t,estim_MC,label='Naive MC',color='blue')\n","plt.fill_between(t, estim_MC-1.96*sigma_MC/nX**.5, estim_MC+1.96*sigma_MC/nX**.5, color='blue', alpha=.1)\n","\n","plt.plot(t,estim_CV,label='Control Variate',color='orange')\n","plt.fill_between(t, estim_CV-1.96*sigma_CV/nX**.5, estim_CV+1.96*sigma_CV/nX**.5, color='orange', alpha=.1)\n","\n","plt.legend()\n","plt.show()"],"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":430},"id":"QiSil_xyNGKX","executionInfo":{"status":"ok","timestamp":1698049633189,"user_tz":-120,"elapsed":574,"user":{"displayName":"Julien Reygner","userId":"06465841509771183939"}},"outputId":"239cde84-6068-486b-9a20-f2d52a60d258"},"execution_count":27,"outputs":[{"output_type":"display_data","data":{"text/plain":["