mirror of
https://github.com/TeachBooks/manual.git
synced 2026-11-12 22:35:51 +00:00
221 lines
7.4 KiB
Text
221 lines
7.4 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Weighted least-squares estimation\n",
|
|
"In ordinary least-squares estimation, we assume that all observations are equally important. In many cases this is not realistic, as observations may be obtained by different measurement systems, or under different circumstances. We want our methodology for least-squares estimation to be able to take this into account.\n",
|
|
"\n",
|
|
"We achieve this goal by introducing a weight matrix in the minimization problem:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\underset{\\mathrm{x}}{\\min} {\\mathrm{(y-Ax)^T W(y-Ax)}}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"In the unweighted least-squares approach, we arrive at the normal equation by pre-multiplying both sides of $\\mathrm{y=Ax}$ with the transpose of the design matrix $\\mathrm{A^T}$:\n",
|
|
"\n",
|
|
"$$ \n",
|
|
"\\mathrm{y=Ax} \\; \\rightarrow \\; \\mathrm{A^T\\; y = A^T\\; A x }\n",
|
|
"$$\n",
|
|
"\n",
|
|
"In the weighted least-squares approach, we now need to add weight matrix $W$ to this pre-multiplication factor, i.e., $ \\mathrm{A^T W}$, to obtain the normal equation:\n",
|
|
"\n",
|
|
"$$ \n",
|
|
"\\mathrm{y=Ax} \\; \\rightarrow \\; \\mathrm{A^T W \\; y = A^TW \\; A x}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"The normal matrix is now defined as $\\mathrm{N=A^T W A}$. From this, assuming that the normal matrix $N$ is invertible (non-singular) we find the weighted least-squares estimate $ \\mathrm{\\hat{x}} $, \n",
|
|
"\n",
|
|
"$$ \n",
|
|
"\\begin{align*}\n",
|
|
"\\mathrm{\\hat{x}} &= \\mathrm{(A^T W A)^{-1} A^T W y} \\\\\n",
|
|
" &= \\arg \\underset{\\mathrm{x}}{\\min} {\\mathrm{(y-Ax)^T W(y-Ax)}}\n",
|
|
"\\end{align*}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"We also find the derived estimate $ \\mathrm{\\hat{y}} $ and $ \\mathrm{\\hat{\\epsilon}} $:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\mathrm{\\hat{y} = A \\hat{x} = A (A^T W A )^{-1} A^T W y}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"$$\n",
|
|
" \\mathrm{\\hat{\\epsilon} = y - \\hat{y}= y - A \\hat{x} = y-A (A^T W A )^{-1} A^T W y = (I- A(A^T W A )^{-1} A^T W) y}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"### Quiz question\n",
|
|
"<iframe src=\"https://tudelft.h5p.com/content/1292046737060674407/embed\" aria-label=\"WLS_1\" width=\"1088\" height=\"637\" frameborder=\"0\" allowfullscreen=\"allowfullscreen\" allow=\"autoplay *; geolocation *; microphone *; camera *; midi *; encrypted-media *\"></iframe><script src=\"https://tudelft.h5p.com/js/h5p-resizer.js\" charset=\"UTF-8\"></script>\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"### Video\n",
|
|
"<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/iJmkkz37EuU\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>\n",
|
|
"\n",
|
|
"### Discussion on the weight matrix\n",
|
|
"The weight matrix $\\mathrm{W}$ expresses the (relative) weights between the observations. It is always a square matrix. The size of the weight matrix depends on the number of observations, $m$. The size of the weight matrix is $m\\times m$. \n",
|
|
"\n",
|
|
"If it is a unit matrix ($\\mathrm{W=I}$), this implies that all observations have equal weight. Note that in this case the equations are equal to the ordinary least-squares solution. \n",
|
|
"\n",
|
|
"If it is a diagonal matrix, with different values on the diagonal, this implies that entries with a higher value correspond to the observations which are considered to be of more importance. If the weight matrix has non-zero elements on the off-diagonal positions, this implies that (some) observations are correlated.\n",
|
|
"\n",
|
|
"### Weighted least-squares estimator: properties\n",
|
|
"\n",
|
|
"Until now, we have looked at the weighted least-squares solution of a single *realization* of the observations, where generally we assume that it is a realization of the *random* observable vector $Y$, since measurements are affected by random errors. As such it follows the the weighted least-squares *estimator* is given by:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\hat{X} = \\mathrm{(A^T W A )^{-1} A^T W} Y\n",
|
|
"$$\n",
|
|
"\n",
|
|
"This estimator has two important properties: it is *linear* and *unbiased*.\n",
|
|
"\n",
|
|
"The linearity property is due to the fact that $\\hat{X}$ is a linear function of the observables $Y$.\n",
|
|
"\n",
|
|
"The unbiased property means that the expectation of $\\hat{X}$ is equal to the true (but unknown) $\\mathrm{x}$. This can be shown as follows:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\mathbb{E}(\\hat{X}) = \\mathrm{(A^T W A )^{-1} A^T W} \\mathbb{E}(Y) = \\mathrm{(A^T W A )^{-1} A^T W Ax = x}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"This a very desirable property. It applies that if we would repeat the measurements many times to obtain a new estimate, the *average of the estimated* values would be equal to the truy values."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Exercise\n",
|
|
"\n",
|
|
"You have a time series of 8 measurements and fit a model assuming a linear trend (constant velocity)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# times of observation [months]\n",
|
|
"t = np.arange(8)\n",
|
|
"\n",
|
|
"# observed heights [m]\n",
|
|
"y = [1.39, 1.26, 1.48, 4.03, 5.89, 5.14, 6.47, 7.64]\n",
|
|
"\n",
|
|
"# number of observations\n",
|
|
"m = len(t)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Fill in the correct $\\mathrm{A}$-matrix."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# design matrix\n",
|
|
"A = np.column_stack((np.ones(m), t))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Define the weight matrix for $\\mathrm{W}=I_m$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Weight matrix for case 1\n",
|
|
"W_1 = ?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Define the weight matrix with the weight of first 3 observations is five times as large as the rest"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Weight matrix for case 2\n",
|
|
"W_2 = ?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Define the weight matrix with the weight of 4th and 5th observation is five times as large as the rest"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Weight matrix for case 3\n",
|
|
"W_3 = ?"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "base",
|
|
"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.9.18"
|
|
},
|
|
"orig_nbformat": 4
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|