{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# AM783: Applied Markov processes | Statistical or Monte Carlo estimation\n",
"\n",
"Hugo Touchette\n",
"\n",
"13 August 2020\n",
"\n",
"Last updated: 11 August 2022\n",
"\n",
"Python 3"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Magic command for vectorised figures\n",
"%config InlineBackend.figure_format = 'svg'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Point estimator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to estimate the mean $E[R]$ of a sample of values known to be distributed from the Rayleigh distribution\n",
"$$\n",
"p(r) = r e^{-r^2/2}.\n",
"$$\n",
"For this, we generate a sample of $L$ values and calculate the estimator\n",
"$$\n",
"\\hat\\mu_L=\\frac{1}{L}\\sum_{i=1}^L X_i.\n",
"$$\n",
"By increasing $L$, the estimator should get closer to the theoretical expectation, which is known to be\n",
"$$\n",
"E[R] = \\int_0^\\infty r\\, p(r)\\, dr = \\sqrt{\\frac{\\pi}{2}}.\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimator value = 1.2525155110475363\n",
"Theoretical value = 1.2533141373155001\n"
]
}
],
"source": [
"L = 10**6\n",
"sample = np.random.rayleigh(1, L) # Using the built-in Rayleigh generator\n",
"mean_est = np.mean(sample)\n",
"print('Estimator value = ', mean_est)\n",
"print('Theoretical value = ', np.sqrt(np.pi/2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Repeating the calculation, we indeed see that estimated values get close to the theoretical expectation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gaussian error bars"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the sample used to obtain the mean estimate, we can also calculate the Gaussian error, as shown in the lecture notes."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimator value = 1.252569432127885 ± 0.0006548888077747249\n",
"Absolute error = 0.0007447051876150468\n"
]
}
],
"source": [
"L = 10**6\n",
"sample = np.random.rayleigh(1, L) # Using the built-in Rayleigh generator\n",
"mean_est = np.mean(sample)\n",
"mean_est_err = np.std(sample)/np.sqrt(L)\n",
"print('Estimator value = ', mean_est, \"±\", mean_est_err)\n",
"print('Absolute error = ', np.abs(mean_est-np.sqrt(np.pi/2)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convergence analysis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An estimation should be repeated for increasing values of $L$ to show that it converges. The code next does this efficiently without recalculating the estimator and its error for each $L$."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"