From e6a7f9704dd7a3b5793057c08bcbf31f53ce1369 Mon Sep 17 00:00:00 2001
From: Kees van Kempen <k.vankempen@student.science.ru.nl>
Date: Tue, 27 Sep 2022 12:21:10 +0200
Subject: [PATCH] 03: Try f_Z(z) = sec^2(pi*z)

---
 Exercise sheet 3/exercise_sheet_03.ipynb | 46 +++++++++++++++++++-----
 1 file changed, 37 insertions(+), 9 deletions(-)

diff --git a/Exercise sheet 3/exercise_sheet_03.ipynb b/Exercise sheet 3/exercise_sheet_03.ipynb
index 3495c43..0b9fa94 100644
--- a/Exercise sheet 3/exercise_sheet_03.ipynb	
+++ b/Exercise sheet 3/exercise_sheet_03.ipynb	
@@ -361,7 +361,7 @@
      "text": [
       "\n",
       "After trial and error, I found that 43402 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
-      "with mean 0.486440 and standard deviation sigma 0.000998. Wolfram Alpha tells me it should be approximately I =\n",
+      "with mean 0.486091 and standard deviation sigma 0.001000. Wolfram Alpha tells me it should be approximately I =\n",
       "0.487595, so I am happy.\n",
       "      \n"
      ]
@@ -390,9 +390,20 @@
     "__(b)__ Choose a random variable $Z$ on $(0,1)$ whose density resembles the integrand of $I$ and which you know how to sample efficiently (by inversion method, acceptance-rejection, or a built-in Python function). Estimate $I$ again using importance sampling, i.e. $I = \\mathbb{E}[X']$ where $X' = g(Z) f_U(Z)/f_Z(Z)$, with an error of at most 0.001. How many samples did you need this time? **(20 pts)**"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "605f8777",
+   "metadata": {},
+   "source": [
+    "After hours of thinking, I tried $$f_Z(x) = \\sec^2(\\pi{}x) = \\frac{1}{\\cos^2(\\pi{}x)}$$ on $(0, 1)$.\n",
+    "It has cumulative density $$F_Z(x) = \\int_0^x\\sec^2(\\pi{}t)dt = \\left[ \\frac{\\tan (\\pi{}t)}{\\pi} \\right]_{t=0}^x = \\frac{\\tan (\\pi{}x)}{\\pi}.$$ This density gives us $$X' := g(Z)f_U(Z)/f_Z(Z) = \\sin(\\pi{}Z(1-Z))\\sin^2(\\pi{}Z)$$ for $Z$ on $(0, 1)$. $Z$ can be sampled using the inversion method as $$F_Z^{-1}(p) = \\frac{\\arctan(\\pi{}p)}{\\pi}.$$\n",
+    "\n",
+    "Unluckily, I found that the function is way too non-constant."
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 8,
    "id": "5793420a",
    "metadata": {
     "deletable": false,
@@ -408,20 +419,37 @@
      "task": false
     }
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "After trial and error, I found that 43402 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
+      "with mean 0.175938 and standard deviation sigma 0.000418. Wolfram Alpha tells me it should be approximately I =\n",
+      "0.487595, so I am happy.\n",
+      "      \n"
+     ]
+    }
+   ],
    "source": [
     "def sample_nice_Z():\n",
     "    '''Sample from the nice distribution Z'''\n",
-    "    # YOUR CODE HERE\n",
-    "    raise NotImplementedError()\n",
+    "    p = rng.random()\n",
+    "    return np.arctan(p*np.pi)/np.pi\n",
     "    \n",
     "def sample_X_prime():\n",
     "    '''Sample from X'.'''\n",
-    "    # YOUR CODE HERE\n",
-    "    raise NotImplementedError()\n",
+    "    z = sample_nice_Z()\n",
+    "    return np.sin(np.pi*z*(1 - z))*np.cos(np.pi*z)**2\n",
     "    \n",
-    "# YOUR CODE HERE\n",
-    "raise NotImplementedError()"
+    "n = 43402\n",
+    "sample_mean, sample_stdev = estimate_expectation(sample_X_prime, n)\n",
+    "print(\"\"\"\n",
+    "After trial and error, I found that {} trials were sufficient to arrive at an 1-sigma error of around 0.001,\n",
+    "with mean {:.6f} and standard deviation sigma {:.6f}. Wolfram Alpha tells me it should be approximately I =\n",
+    "{}, so I am happy.\n",
+    "      \"\"\".format(n, sample_mean, sample_stdev, 0.487595))"
    ]
   },
   {