diff --git a/doc/pub/Regression/ipynb/Regression.ipynb b/doc/pub/Regression/ipynb/Regression.ipynb
index 4bc47575b..65a7314ca 100644
--- a/doc/pub/Regression/ipynb/Regression.ipynb
+++ b/doc/pub/Regression/ipynb/Regression.ipynb
@@ -381,7 +381,9 @@
{
"cell_type": "code",
"execution_count": 1,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"%matplotlib inline\n",
@@ -886,7 +888,9 @@
{
"cell_type": "code",
"execution_count": 2,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# matrix inversion to find beta\n",
@@ -905,7 +909,9 @@
{
"cell_type": "code",
"execution_count": 3,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fit = np.linalg.lstsq(X, Energies, rcond =None)[0]\n",
@@ -922,7 +928,9 @@
{
"cell_type": "code",
"execution_count": 4,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"Masses['Eapprox'] = ytilde\n",
@@ -952,7 +960,9 @@
{
"cell_type": "code",
"execution_count": 5,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def R2(y_data, y_model):\n",
@@ -969,7 +979,9 @@
{
"cell_type": "code",
"execution_count": 6,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"print(R2(Energies,ytilde))"
@@ -985,7 +997,9 @@
{
"cell_type": "code",
"execution_count": 7,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def MSE(y_data,y_model):\n",
@@ -1005,7 +1019,9 @@
{
"cell_type": "code",
"execution_count": 8,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def RelativeError(y_data,y_model):\n",
@@ -1382,7 +1398,9 @@
{
"cell_type": "code",
"execution_count": 9,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Common imports\n",
@@ -1501,7 +1519,9 @@
{
"cell_type": "code",
"execution_count": 10,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import os\n",
@@ -1618,7 +1638,9 @@
{
"cell_type": "code",
"execution_count": 11,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -1638,7 +1660,9 @@
{
"cell_type": "code",
"execution_count": 12,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from sklearn.datasets import load_boston\n",
@@ -1660,7 +1684,9 @@
{
"cell_type": "code",
"execution_count": 13,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"boston = pd.DataFrame(boston_dataset.data, columns=boston_dataset.feature_names)\n",
@@ -1678,7 +1704,9 @@
{
"cell_type": "code",
"execution_count": 14,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# check for missing values in all the columns\n",
@@ -1695,7 +1723,9 @@
{
"cell_type": "code",
"execution_count": 15,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# set the size of the figure\n",
@@ -1716,7 +1746,9 @@
{
"cell_type": "code",
"execution_count": 16,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# compute the pair wise correlation for all columns \n",
@@ -1736,7 +1768,9 @@
{
"cell_type": "code",
"execution_count": 17,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"plt.figure(figsize=(20, 5))\n",
@@ -1764,7 +1798,9 @@
{
"cell_type": "code",
"execution_count": 18,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X = pd.DataFrame(np.c_[boston['LSTAT'], boston['RM']], columns = ['LSTAT','RM'])\n",
@@ -1781,7 +1817,9 @@
{
"cell_type": "code",
"execution_count": 19,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from sklearn.model_selection import train_test_split\n",
@@ -1805,7 +1843,9 @@
{
"cell_type": "code",
"execution_count": 20,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from sklearn.linear_model import LinearRegression\n",
@@ -1844,7 +1884,9 @@
{
"cell_type": "code",
"execution_count": 21,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# plotting the y_test vs y_pred\n",
@@ -1929,7 +1971,9 @@
{
"cell_type": "code",
"execution_count": 22,
- "metadata": {},
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Common imports\n",
@@ -2045,7 +2089,7 @@
"however not the be case in general and a standard matrix inversion\n",
"algorithm based on say LU, QR or Cholesky decomposition may lead to singularities. We will see examples of this below.\n",
"\n",
- "There is however a way to partially circumvent this problem and also gain some insight about the ordinary least squares approach. \n",
+ "There is however a way to partially circumvent this problem and also gain some insights about the ordinary least squares approach, and later shrinkage methods like Ridge and Lasso regressions. \n",
"\n",
"This is given by the **Singular Value Decomposition** algorithm, perhaps\n",
"the most powerful linear algebra algorithm. Let us look at a\n",
@@ -2299,32 +2343,9 @@
"with eigenvalues $\\sigma_1=2$ and $\\sigma_2=0$. \n",
"The SVD exits always! \n",
"\n",
- "\n",
- "## Another Example\n",
- "\n",
- "Consider the following matrix which can be SVD decomposed as"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\boldsymbol{X} = \\frac{1}{15}\\begin{bmatrix} 14 & 2\\\\ 4 & 22\\\\ 16 & 13\\end{bmatrix}=\\frac{1}{3}\\begin{bmatrix} 1& 2 & 2 \\\\ 2& -1 & 1\\\\ 2 & 1& -2\\end{bmatrix} \\begin{bmatrix} 2& 0 \\\\ 0& 1\\\\ 0 & 0\\end{bmatrix}\\frac{1}{5}\\begin{bmatrix} 3& 4 \\\\ 4& -3\\end{bmatrix}=\\boldsymbol{U}\\boldsymbol{\\Sigma}\\boldsymbol{V}^T.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "This is a $3\\times 2$ matrix which is decomposed in terms of a\n",
- "$3\\times 3$ matrix $\\boldsymbol{U}$, and a $2\\times 2$ matrix $\\boldsymbol{V}$. It is easy to see\n",
- "that $\\boldsymbol{U}$ and $\\boldsymbol{V}$ are orthogonal (how?). \n",
- "\n",
- "And the SVD\n",
+ "The SVD\n",
"decomposition (singular values) gives eigenvalues \n",
- "$\\sigma_i\\geq\\sigma_{i+1}$ for all $i$ and for dimensions larger than $i=2$, the\n",
+ "$\\sigma_i\\geq\\sigma_{i+1}$ for all $i$ and for dimensions larger than $i=p$, the\n",
"eigenvalues (singular values) are zero.\n",
"\n",
"In the general case, where our design matrix $\\boldsymbol{X}$ has dimension\n",
@@ -2359,6 +2380,66 @@
"The $n=p$ case is obvious, we retain the full SVD. \n",
"In general the economy-size SVD leads to less FLOPS and still conserving the desired accuracy.\n",
"\n",
+ "## Codes for the SVD"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "# SVD inversion\n",
+ "def SVDinv(A):\n",
+ " ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).\n",
+ " SVD is numerically more stable than the inversion algorithms provided by\n",
+ " numpy and scipy.linalg at the cost of being slower.\n",
+ " '''\n",
+ " U, s, VT = np.linalg.svd(A)\n",
+ "# print('test U')\n",
+ "# print( (np.transpose(U) @ U - U @np.transpose(U)))\n",
+ "# print('test VT')\n",
+ "# print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))\n",
+ " print(U)\n",
+ " print(s)\n",
+ " print(VT)\n",
+ "\n",
+ " D = np.zeros((len(U),len(VT)))\n",
+ " for i in range(0,len(VT)):\n",
+ " D[i,i]=s[i]\n",
+ " UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)\n",
+ " return np.matmul(V,np.matmul(invD,UT))\n",
+ "\n",
+ "\n",
+ "X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])\n",
+ "print(X)\n",
+ "A = np.transpose(X) @ X\n",
+ "print(A)\n",
+ "# Brute force inversion of super-collinear matrix\n",
+ "#B = np.linalg.inv(A)\n",
+ "#print(B)\n",
+ "C = SVDinv(A)\n",
+ "print(C)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The matrix $\\boldsymbol{X}$ has columns that are linearly dependent. The first\n",
+ "column is the row-wise sum of the other two columns. The rank of a\n",
+ "matrix (the column rank) is the dimension of space spanned by the\n",
+ "column vectors. The rank of the matrix is the number of linearly\n",
+ "independent columns, in this case just $2$. We see this from the\n",
+ "singular values when running the above code. Running the standard\n",
+ "inversion algorithm for matrix inversion with $\\boldsymbol{X}^T\\boldsymbol{X}$ results\n",
+ "in the program terminating due to a singular matrix.\n",
+ "\n",
+ "\n",
+ "\n",
"## Mathematical Properties\n",
"\n",
"There are several interesting mathematical properties which will be\n",
@@ -2790,86 +2871,6 @@
"For more discussions of Ridge and Lasso regression, [Wessel van Wieringen's](https://arxiv.org/abs/1509.09169) article is highly recommended.\n",
"Similarly, [Mehta et al's article](https://arxiv.org/abs/1803.08823) is also recommended.\n",
"\n",
- "## Codes for the SVD"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[[ 1. -1. 2.]\n",
- " [ 1. 0. 1.]\n",
- " [ 1. 2. -1.]\n",
- " [ 1. 1. 0.]]\n",
- "[[ 4. 2. 2.]\n",
- " [ 2. 6. -4.]\n",
- " [ 2. -4. 6.]]\n",
- "[[-9.57425734e-17 8.16496581e-01 -5.77350269e-01]\n",
- " [-7.07106781e-01 4.08248290e-01 5.77350269e-01]\n",
- " [ 7.07106781e-01 4.08248290e-01 5.77350269e-01]]\n",
- "[1.00000000e+01 6.00000000e+00 9.10898112e-32]\n",
- "[[ 3.33066907e-17 -7.07106781e-01 7.07106781e-01]\n",
- " [ 8.16496581e-01 4.08248290e-01 4.08248290e-01]\n",
- " [ 5.77350269e-01 -5.77350269e-01 -5.77350269e-01]]\n",
- "[[-3.65939208e+30 3.65939208e+30 3.65939208e+30]\n",
- " [ 3.65939208e+30 -3.65939208e+30 -3.65939208e+30]\n",
- " [ 3.65939208e+30 -3.65939208e+30 -3.65939208e+30]]\n"
- ]
- }
- ],
- "source": [
- "import numpy as np\n",
- "# SVD inversion\n",
- "def SVDinv(A):\n",
- " ''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).\n",
- " SVD is numerically more stable than the inversion algorithms provided by\n",
- " numpy and scipy.linalg at the cost of being slower.\n",
- " '''\n",
- " U, s, VT = np.linalg.svd(A)\n",
- "# print('test U')\n",
- "# print( (np.transpose(U) @ U - U @np.transpose(U)))\n",
- "# print('test VT')\n",
- "# print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))\n",
- " print(U)\n",
- " print(s)\n",
- " print(VT)\n",
- "\n",
- " D = np.zeros((len(U),len(VT)))\n",
- " for i in range(0,len(VT)):\n",
- " D[i,i]=s[i]\n",
- " UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)\n",
- " return np.matmul(V,np.matmul(invD,UT))\n",
- "\n",
- "\n",
- "X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])\n",
- "print(X)\n",
- "A = np.transpose(X) @ X\n",
- "print(A)\n",
- "# Brute force inversion of super-collinear matrix\n",
- "#B = np.linalg.inv(A)\n",
- "#print(B)\n",
- "C = SVDinv(A)\n",
- "print(C)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The matrix $\\boldsymbol{X}$ has columns that are linearly dependent. The first\n",
- "column is the row-wise sum of the other two columns. The rank of a\n",
- "matrix (the column rank) is the dimension of space spanned by the\n",
- "column vectors. The rank of the matrix is the number of linearly\n",
- "independent columns, in this case just $2$. We see this from the\n",
- "singular values when running the above code. Running the standard\n",
- "inversion algorithm for matrix inversion with $\\boldsymbol{X}^T\\boldsymbol{X}$ results\n",
- "in the program terminating due to a singular matrix.\n",
- "\n",
"\n",
"\n",
"## A better understanding of regularization\n",
@@ -3176,20 +3177,11 @@
},
{
"cell_type": "code",
- "execution_count": 2,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "0.06921280734901214\n",
- "4.282897881229527\n",
- "[[ 1.16689268 3.45320489]\n",
- " [ 3.45320489 11.11096825]]\n"
- ]
- }
- ],
+ "execution_count": 24,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"# Importing various packages\n",
"import numpy as np\n",
@@ -3218,20 +3210,11 @@
},
{
"cell_type": "code",
- "execution_count": 3,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "0.07629059232706191\n",
- "1.3689323216991471\n",
- "[[1. 0.63967543]\n",
- " [0.63967543 1. ]]\n"
- ]
- }
- ],
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"import numpy as np\n",
"n = 100\n",
@@ -3273,40 +3256,11 @@
},
{
"cell_type": "code",
- "execution_count": 4,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[[-0.55783318 -2.36855737]\n",
- " [-0.48038212 -1.12123656]\n",
- " [-0.90619932 -1.43972046]\n",
- " [ 0.0509597 0.25649197]\n",
- " [-0.00589746 1.0155279 ]\n",
- " [ 0.42557278 1.62052588]\n",
- " [ 0.76914274 -1.09850401]\n",
- " [-0.13134882 0.66162568]\n",
- " [ 0.85273632 2.51549691]\n",
- " [-0.01675064 -0.04164992]]\n",
- " 0 1\n",
- "0 -0.557833 -2.368557\n",
- "1 -0.480382 -1.121237\n",
- "2 -0.906199 -1.439720\n",
- "3 0.050960 0.256492\n",
- "4 -0.005897 1.015528\n",
- "5 0.425573 1.620526\n",
- "6 0.769143 -1.098504\n",
- "7 -0.131349 0.661626\n",
- "8 0.852736 2.515497\n",
- "9 -0.016751 -0.041650\n",
- " 0 1\n",
- "0 1.000000 0.657956\n",
- "1 0.657956 1.000000\n"
- ]
- }
- ],
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
@@ -3334,49 +3288,11 @@
},
{
"cell_type": "code",
- "execution_count": 5,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- " 0 1 2 3 4 5 6 7 \\\n",
- "0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n",
- "1 0.0 0.098274 0.088019 0.102816 0.095130 0.088038 0.095734 0.089350 \n",
- "2 0.0 0.088019 0.080153 0.091357 0.085370 0.079873 0.085008 0.079956 \n",
- "3 0.0 0.102816 0.091357 0.113039 0.104455 0.096480 0.108457 0.101346 \n",
- "4 0.0 0.095130 0.085370 0.104455 0.097137 0.090339 0.100399 0.094305 \n",
- "5 0.0 0.088038 0.079873 0.096480 0.090339 0.084640 0.092879 0.087723 \n",
- "6 0.0 0.095734 0.085008 0.108457 0.100399 0.092879 0.106219 0.099527 \n",
- "7 0.0 0.089350 0.079956 0.101346 0.094305 0.087723 0.099527 0.093665 \n",
- "8 0.0 0.083518 0.075349 0.094820 0.088708 0.082989 0.093363 0.088259 \n",
- "9 0.0 0.078178 0.071145 0.088809 0.083553 0.078629 0.087665 0.083254 \n",
- "10 0.0 0.087401 0.077865 0.101048 0.093860 0.087128 0.100476 0.094478 \n",
- "11 0.0 0.081937 0.073490 0.094954 0.088611 0.082657 0.094714 0.089415 \n",
- "12 0.0 0.076958 0.069504 0.089383 0.083807 0.078560 0.089433 0.084768 \n",
- "13 0.0 0.072416 0.065869 0.084280 0.079402 0.074802 0.084581 0.080492 \n",
- "14 0.0 0.068264 0.062550 0.079595 0.075355 0.071347 0.080112 0.076548 \n",
- "\n",
- " 8 9 10 11 12 13 14 \n",
- "0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 \n",
- "1 0.083518 0.078178 0.087401 0.081937 0.076958 0.072416 0.068264 \n",
- "2 0.075349 0.071145 0.077865 0.073490 0.069504 0.065869 0.062550 \n",
- "3 0.094820 0.088809 0.101048 0.094954 0.089383 0.084280 0.079595 \n",
- "4 0.088708 0.083553 0.093860 0.088611 0.083807 0.079402 0.075355 \n",
- "5 0.082989 0.078629 0.087128 0.082657 0.078560 0.074802 0.071347 \n",
- "6 0.093363 0.087665 0.100476 0.094714 0.089433 0.084581 0.080112 \n",
- "7 0.088259 0.083254 0.094478 0.089415 0.084768 0.080492 0.076548 \n",
- "8 0.083545 0.079175 0.088938 0.084512 0.080443 0.076693 0.073228 \n",
- "9 0.079175 0.075391 0.083801 0.079956 0.076416 0.073148 0.070125 \n",
- "10 0.088938 0.083801 0.096178 0.090983 0.086211 0.081816 0.077756 \n",
- "11 0.084512 0.079956 0.090983 0.086383 0.082149 0.078242 0.074627 \n",
- "12 0.080443 0.076416 0.086211 0.082149 0.078403 0.074940 0.071730 \n",
- "13 0.076693 0.073148 0.081816 0.078242 0.074940 0.071881 0.069039 \n",
- "14 0.073228 0.070125 0.077756 0.074627 0.071730 0.069039 0.066534 \n"
- ]
- }
- ],
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
"source": [
"# Common imports\n",
"import numpy as np\n",
@@ -3484,7 +3400,7 @@
"metadata": {},
"source": [
"$$\n",
- "\\mathbb{E}[\\boldsymbol{X}\\boldsymbol{X}^T] = \\frac{1}{n}\\boldsymbol{X}\\boldsymbol{X}^T=\\begin{bmatrix}\n",
+ "\\mathbb{E}[\\boldsymbol{X}^T\\boldsymbol{X}] = \\frac{1}{n}\\boldsymbol{X}^T\\boldsymbol{X}=\\begin{bmatrix}\n",
"x_{00}^2+x_{01}^2 & x_{00}x_{10}+x_{01}x_{11}\\\\\n",
"x_{10}x_{00}+x_{11}x_{01} & x_{10}^2+x_{11}^2\\\\\n",
"\\end{bmatrix},\n",
@@ -3612,1265 +3528,11 @@
"\n",
"\n",
"\n",
- "## Statistics\n",
- "The *probability distribution function (PDF)* is a function\n",
- "$p(x)$ on the domain which, in the discrete case, gives us the\n",
- "probability or relative frequency with which these values of $X$ occur:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "p(x) = \\mathrm{prob}(X=x)\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "In the continuous case, the PDF does not directly depict the\n",
- "actual probability. Instead we define the probability for the\n",
- "stochastic variable to assume any value on an infinitesimal interval\n",
- "around $x$ to be $p(x)dx$. The continuous function $p(x)$ then gives us\n",
- "the *density* of the probability rather than the probability\n",
- "itself. The probability for a stochastic variable to assume any value\n",
- "on a non-infinitesimal interval $[a,\\,b]$ is then just the integral:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{prob}(a\\leq X\\leq b) = \\int_a^b p(x)dx\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Qualitatively speaking, a stochastic variable represents the values of\n",
- "numbers chosen as if by chance from some specified PDF so that the\n",
- "selection of a large set of these numbers reproduces this PDF.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, moments\n",
- "A particularly useful class of special expectation values are the\n",
- "*moments*. The $n$-th moment of the PDF $p$ is defined as\n",
- "follows:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\langle x^n\\rangle \\equiv \\int\\! x^n p(x)\\,dx\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The zero-th moment $\\langle 1\\rangle$ is just the normalization condition of\n",
- "$p$. The first moment, $\\langle x\\rangle$, is called the *mean* of $p$\n",
- "and often denoted by the letter $\\mu$:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\langle x\\rangle = \\mu \\equiv \\int\\! x p(x)\\,dx\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, central moments\n",
- "A special version of the moments is the set of *central moments*,\n",
- "the n-th central moment defined as:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\langle (x-\\langle x \\rangle )^n\\rangle \\equiv \\int\\! (x-\\langle x\\rangle)^n p(x)\\,dx\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The zero-th and first central moments are both trivial, equal $1$ and\n",
- "$0$, respectively. But the second central moment, known as the\n",
- "*variance* of $p$, is of particular interest. For the stochastic\n",
- "variable $X$, the variance is denoted as $\\sigma^2_X$ or $\\mathrm{var}(X)$:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\sigma^2_X\\ \\ =\\ \\ \\mathrm{var}(X) = \\langle (x-\\langle x\\rangle)^2\\rangle =\n",
- "\\int\\! (x-\\langle x\\rangle)^2 p(x)\\,dx\n",
- "\\label{_auto2} \\tag{2}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- " = \\int\\! \\left(x^2 - 2 x \\langle x\\rangle^{2} +\n",
- " \\langle x\\rangle^2\\right)p(x)\\,dx\n",
- "\\label{_auto3} \\tag{3}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- " = \\langle x^2\\rangle - 2 \\langle x\\rangle\\langle x\\rangle + \\langle x\\rangle^2\n",
- "\\label{_auto4} \\tag{4}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- " = \\langle x^2\\rangle - \\langle x\\rangle^2\n",
- "\\label{_auto5} \\tag{5}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The square root of the variance, $\\sigma =\\sqrt{\\langle (x-\\langle x\\rangle)^2\\rangle}$ is called the *standard deviation* of $p$. It is clearly just the RMS (root-mean-square)\n",
- "value of the deviation of the PDF from its mean value, interpreted\n",
- "qualitatively as the *spread* of $p$ around its mean.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, covariance\n",
- "Another important quantity is the so called covariance, a variant of\n",
- "the above defined variance. Consider again the set $\\{X_i\\}$ of $n$\n",
- "stochastic variables (not necessarily uncorrelated) with the\n",
- "multivariate PDF $P(x_1,\\dots,x_n)$. The *covariance* of two\n",
- "of the stochastic variables, $X_i$ and $X_j$, is defined as follows:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{cov}(X_i,\\,X_j) \\equiv \\langle (x_i-\\langle x_i\\rangle)(x_j-\\langle x_j\\rangle)\\rangle\n",
- "\\nonumber\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\n",
- "\\int\\!\\cdots\\!\\int\\!(x_i-\\langle x_i \\rangle)(x_j-\\langle x_j \\rangle)\\,\n",
- "P(x_1,\\dots,x_n)\\,dx_1\\dots dx_n\n",
- "\\label{eq:def_covariance} \\tag{6}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "with"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\langle x_i\\rangle =\n",
- "\\int\\!\\cdots\\!\\int\\!x_i\\,P(x_1,\\dots,x_n)\\,dx_1\\dots dx_n\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, more covariance\n",
- "If we consider the above covariance as a matrix $C_{ij}=\\mathrm{cov}(X_i,\\,X_j)$, then the diagonal elements are just the familiar\n",
- "variances, $C_{ii} = \\mathrm{cov}(X_i,\\,X_i) = \\mathrm{var}(X_i)$. It turns out that\n",
- "all the off-diagonal elements are zero if the stochastic variables are\n",
- "uncorrelated. This is easy to show, keeping in mind the linearity of\n",
- "the expectation value. Consider the stochastic variables $X_i$ and\n",
- "$X_j$, ($i\\neq j$):"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\mathrm{cov}(X_i,\\,X_j) = \\langle(x_i-\\langle x_i\\rangle)(x_j-\\langle x_j\\rangle)\\rangle\n",
- "\\label{_auto6} \\tag{7}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\langle x_i x_j - x_i\\langle x_j\\rangle - \\langle x_i\\rangle x_j + \\langle x_i\\rangle\\langle x_j\\rangle\\rangle \n",
- "\\label{_auto7} \\tag{8}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\langle x_i x_j\\rangle - \\langle x_i\\langle x_j\\rangle\\rangle - \\langle \\langle x_i\\rangle x_j\\rangle +\n",
- "\\langle \\langle x_i\\rangle\\langle x_j\\rangle\\rangle\n",
- "\\label{_auto8} \\tag{9}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\langle x_i x_j\\rangle - \\langle x_i\\rangle\\langle x_j\\rangle - \\langle x_i\\rangle\\langle x_j\\rangle +\n",
- "\\langle x_i\\rangle\\langle x_j\\rangle\n",
- "\\label{_auto9} \\tag{10}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\langle x_i x_j\\rangle - \\langle x_i\\rangle\\langle x_j\\rangle\n",
- "\\label{_auto10} \\tag{11}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Covariance example\n",
- "\n",
- "Suppose we have defined three vectors $\\boldsymbol{x}, \\boldsymbol{y}, \\boldsymbol{z}$ with\n",
- "$n$ elements each. The covariance matrix is defined as"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\boldsymbol{\\Sigma} = \\begin{bmatrix} \\sigma_{xx} & \\sigma_{xy} & \\sigma_{xz} \\\\\n",
- " \\sigma_{yx} & \\sigma_{yy} & \\sigma_{yz} \\\\\n",
- " \\sigma_{zx} & \\sigma_{zy} & \\sigma_{zz}\n",
- " \\end{bmatrix},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where for example"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\sigma_{xy} =\\frac{1}{n} \\sum_{i=0}^{n-1}(x_i- \\overline{x})(y_i- \\overline{y}).\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The Numpy function **np.cov** calculates the covariance elements using\n",
- "the factor $1/(n-1)$ instead of $1/n$ since it assumes we do not have\n",
- "the exact mean valu\\ es.\n",
- "\n",
- "The following simple function uses the **np.vstack** function which\n",
- "takes each vector of dimension $1\\times n$ and produces a $3\\times n$\n",
- "matrix $\\boldsymbol{W}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\boldsymbol{W} = \\begin{bmatrix} x_0 & y_0 & z_0 \\\\\n",
- " x_1 & y_1 & z_1 \\\\\n",
- " x_2 & y_2 & z_2 \\\\\n",
- " \\dots & \\dots & \\dots \\\\\n",
- " x_{n-2} & y_{n-2} & z_{n-2} \\\\\n",
- " x_{n-1} & y_{n-1} & z_{n-1}\n",
- " \\end{bmatrix},\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which in turn is converted into into the $3\\times 3$ covariance matrix\n",
- "$\\boldsymbol{\\Sigma}$ via the Numpy function **np.cov()**. We note that we can\n",
- "also calculate the mean value of each set of samples $\\boldsymbol{x}$ etc\n",
- "using the Numpy function **np.mean(x)**. We can also extract the\n",
- "eigenvalues of the covariance matrix through the **np.linalg.eig()**\n",
- "function.\n",
- "\n",
- "\n",
- "## Covariance in numpy"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 28,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Importing various packages\n",
- "import numpy as np\n",
- "\n",
- "n = 100\n",
- "x = np.random.normal(size=n)\n",
- "print(np.mean(x))\n",
- "y = 4+3*x+np.random.normal(size=n)\n",
- "print(np.mean(y))\n",
- "z = x**3+np.random.normal(size=n)\n",
- "print(np.mean(z))\n",
- "W = np.vstack((x, y, z))\n",
- "Sigma = np.cov(W)\n",
- "print(Sigma)\n",
- "Eigvals, Eigvecs = np.linalg.eig(Sigma)\n",
- "print(Eigvals)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 29,
- "metadata": {},
- "outputs": [],
- "source": [
- "import numpy as np\n",
- "import matplotlib.pyplot as plt\n",
- "from scipy import sparse\n",
- "eye = np.eye(4)\n",
- "print(eye)\n",
- "sparse_mtx = sparse.csr_matrix(eye)\n",
- "print(sparse_mtx)\n",
- "x = np.linspace(-10,10,100)\n",
- "y = np.sin(x)\n",
- "plt.plot(x,y,marker='x')\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, independent variables\n",
- "If $X_i$ and $X_j$ are independent, we get \n",
- "$\\langle x_i x_j\\rangle =\\langle x_i\\rangle\\langle x_j\\rangle$, resulting in $\\mathrm{cov}(X_i, X_j) = 0\\ \\ (i\\neq j)$.\n",
- "\n",
- "Also useful for us is the covariance of linear combinations of\n",
- "stochastic variables. Let $\\{X_i\\}$ and $\\{Y_i\\}$ be two sets of\n",
- "stochastic variables. Let also $\\{a_i\\}$ and $\\{b_i\\}$ be two sets of\n",
- "scalars. Consider the linear combination:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "U = \\sum_i a_i X_i \\qquad V = \\sum_j b_j Y_j\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "By the linearity of the expectation value"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{cov}(U, V) = \\sum_{i,j}a_i b_j \\mathrm{cov}(X_i, Y_j)\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, more variance\n",
- "Now, since the variance is just $\\mathrm{var}(X_i) = \\mathrm{cov}(X_i, X_i)$, we get\n",
- "the variance of the linear combination $U = \\sum_i a_i X_i$:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\mathrm{var}(U) = \\sum_{i,j}a_i a_j \\mathrm{cov}(X_i, X_j)\n",
- "\\label{eq:variance_linear_combination} \\tag{12}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "And in the special case when the stochastic variables are\n",
- "uncorrelated, the off-diagonal elements of the covariance are as we\n",
- "know zero, resulting in:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "1\n",
- "2\n",
- "8\n",
- " \n",
- "<\n",
- "<\n",
- "<\n",
- "!\n",
- "!\n",
- "M\n",
- "A\n",
- "T\n",
- "H\n",
- "_\n",
- "B\n",
- "L\n",
- "O\n",
- "C\n",
- "K"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{var}(\\sum_i a_i X_i) = \\sum_i a_i^2 \\mathrm{var}(X_i)\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which will become very useful in our study of the error in the mean\n",
- "value of a set of measurements.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics and stochastic processes\n",
- "A *stochastic process* is a process that produces sequentially a\n",
- "chain of values:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\{x_1, x_2,\\dots\\,x_k,\\dots\\}.\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "We will call these\n",
- "values our *measurements* and the entire set as our measured\n",
- "*sample*. The action of measuring all the elements of a sample\n",
- "we will call a stochastic *experiment* since, operationally,\n",
- "they are often associated with results of empirical observation of\n",
- "some physical or mathematical phenomena; precisely an experiment. We\n",
- "assume that these values are distributed according to some \n",
- "PDF $p_X^{\\phantom X}(x)$, where $X$ is just the formal symbol for the\n",
- "stochastic variable whose PDF is $p_X^{\\phantom X}(x)$. Instead of\n",
- "trying to determine the full distribution $p$ we are often only\n",
- "interested in finding the few lowest moments, like the mean\n",
- "$\\mu_X^{\\phantom X}$ and the variance $\\sigma_X^{\\phantom X}$.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Statistics and sample variables\n",
- "In practical situations a sample is always of finite size. Let that\n",
- "size be $n$. The expectation value of a sample, the *sample mean*, is then defined as follows:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\bar{x}_n \\equiv \\frac{1}{n}\\sum_{k=1}^n x_k\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The *sample variance* is:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{var}(x) \\equiv \\frac{1}{n}\\sum_{k=1}^n (x_k - \\bar{x}_n)^2\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "its square root being the *standard deviation of the sample*. The\n",
- "*sample covariance* is:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{cov}(x)\\equiv\\frac{1}{n}\\sum_{kl}(x_k - \\bar{x}_n)(x_l - \\bar{x}_n)\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, sample variance and covariance\n",
- "Note that the sample variance is the sample covariance without the\n",
- "cross terms. In a similar manner as the covariance in Eq. ([6](#eq:def_covariance)) is a measure of the correlation between\n",
- "two stochastic variables, the above defined sample covariance is a\n",
- "measure of the sequential correlation between succeeding measurements\n",
- "of a sample.\n",
- "\n",
- "These quantities, being known experimental values, differ\n",
- "significantly from and must not be confused with the similarly named\n",
- "quantities for stochastic variables, mean $\\mu_X$, variance $\\mathrm{var}(X)$\n",
- "and covariance $\\mathrm{cov}(X,Y)$.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, law of large numbers\n",
- "The law of large numbers\n",
- "states that as the size of our sample grows to infinity, the sample\n",
- "mean approaches the true mean $\\mu_X^{\\phantom X}$ of the chosen PDF:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\lim_{n\\to\\infty}\\bar{x}_n = \\mu_X^{\\phantom X}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The sample mean $\\bar{x}_n$ works therefore as an estimate of the true\n",
- "mean $\\mu_X^{\\phantom X}$.\n",
- "\n",
- "What we need to find out is how good an approximation $\\bar{x}_n$ is to\n",
- "$\\mu_X^{\\phantom X}$. In any stochastic measurement, an estimated\n",
- "mean is of no use to us without a measure of its error. A quantity\n",
- "that tells us how well we can reproduce it in another experiment. We\n",
- "are therefore interested in the PDF of the sample mean itself. Its\n",
- "standard deviation will be a measure of the spread of sample means,\n",
- "and we will simply call it the *error* of the sample mean, or\n",
- "just sample error, and denote it by $\\mathrm{err}_X^{\\phantom X}$. In\n",
- "practice, we will only be able to produce an *estimate* of the\n",
- "sample error since the exact value would require the knowledge of the\n",
- "true PDFs behind, which we usually do not have.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, more on sample error\n",
- "Let us first take a look at what happens to the sample error as the\n",
- "size of the sample grows. In a sample, each of the measurements $x_i$\n",
- "can be associated with its own stochastic variable $X_i$. The\n",
- "stochastic variable $\\overline X_n$ for the sample mean $\\bar{x}_n$ is\n",
- "then just a linear combination, already familiar to us:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\overline X_n = \\frac{1}{n}\\sum_{i=1}^n X_i\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "All the coefficients are just equal $1/n$. The PDF of $\\overline X_n$,\n",
- "denoted by $p_{\\overline X_n}(x)$ is the desired PDF of the sample\n",
- "means.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics\n",
- "The probability density of obtaining a sample mean $\\bar x_n$\n",
- "is the product of probabilities of obtaining arbitrary values $x_1,\n",
- "x_2,\\dots,x_n$ with the constraint that the mean of the set $\\{x_i\\}$\n",
- "is $\\bar x_n$:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "p_{\\overline X_n}(x) = \\int p_X^{\\phantom X}(x_1)\\cdots\n",
- "\\int p_X^{\\phantom X}(x_n)\\ \n",
- "\\delta\\!\\left(x - \\frac{x_1+x_2+\\dots+x_n}{n}\\right)dx_n \\cdots dx_1\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "And in particular we are interested in its variance $\\mathrm{var}(\\overline X_n)$.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, central limit theorem\n",
- "It is generally not possible to express $p_{\\overline X_n}(x)$ in a\n",
- "closed form given an arbitrary PDF $p_X^{\\phantom X}$ and a number\n",
- "$n$. But for the limit $n\\to\\infty$ it is possible to make an\n",
- "approximation. The very important result is called *the central limit theorem*. It tells us that as $n$ goes to infinity,\n",
- "$p_{\\overline X_n}(x)$ approaches a Gaussian distribution whose mean\n",
- "and variance equal the true mean and variance, $\\mu_{X}^{\\phantom X}$\n",
- "and $\\sigma_{X}^{2}$, respectively:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\lim_{n\\to\\infty} p_{\\overline X_n}(x) =\n",
- "\\left(\\frac{n}{2\\pi\\mathrm{var}(X)}\\right)^{1/2}\n",
- "e^{-\\frac{n(x-\\bar x_n)^2}{2\\mathrm{var}(X)}}\n",
- "\\label{eq:central_limit_gaussian} \\tag{13}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, more technicalities\n",
- "The desired variance\n",
- "$\\mathrm{var}(\\overline X_n)$, i.e. the sample error squared\n",
- "$\\mathrm{err}_X^2$, is given by:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\mathrm{err}_X^2 = \\mathrm{var}(\\overline X_n) = \\frac{1}{n^2}\n",
- "\\sum_{ij} \\mathrm{cov}(X_i, X_j)\n",
- "\\label{eq:error_exact} \\tag{14}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "We see now that in order to calculate the exact error of the sample\n",
- "with the above expression, we would need the true means\n",
- "$\\mu_{X_i}^{\\phantom X}$ of the stochastic variables $X_i$. To\n",
- "calculate these requires that we know the true multivariate PDF of all\n",
- "the $X_i$. But this PDF is unknown to us, we have only got the measurements of\n",
- "one sample. The best we can do is to let the sample itself be an\n",
- "estimate of the PDF of each of the $X_i$, estimating all properties of\n",
- "$X_i$ through the measurements of the sample.\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "## Statistics\n",
- "Our estimate of $\\mu_{X_i}^{\\phantom X}$ is then the sample mean $\\bar x$\n",
- "itself, in accordance with the the central limit theorem:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mu_{X_i}^{\\phantom X} = \\langle x_i\\rangle \\approx \\frac{1}{n}\\sum_{k=1}^n x_k = \\bar x\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Using $\\bar x$ in place of $\\mu_{X_i}^{\\phantom X}$ we can give an\n",
- "*estimate* of the covariance in Eq. ([14](#eq:error_exact))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{cov}(X_i, X_j) = \\langle (x_i-\\langle x_i\\rangle)(x_j-\\langle x_j\\rangle)\\rangle\n",
- "\\approx\\langle (x_i - \\bar x)(x_j - \\bar{x})\\rangle,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "resulting in"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\frac{1}{n} \\sum_{l}^n \\left(\\frac{1}{n}\\sum_{k}^n (x_k -\\bar x_n)(x_l - \\bar x_n)\\right)=\\frac{1}{n}\\frac{1}{n} \\sum_{kl} (x_k -\\bar x_n)(x_l - \\bar x_n)=\\frac{1}{n}\\mathrm{cov}(x)\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics and sample variance\n",
- "By the same procedure we can use the sample variance as an\n",
- "estimate of the variance of any of the stochastic variables $X_i$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{var}(X_i)=\\langle x_i - \\langle x_i\\rangle\\rangle \\approx \\langle x_i - \\bar x_n\\rangle\\nonumber,\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which is approximated as"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\mathrm{var}(X_i)\\approx \\frac{1}{n}\\sum_{k=1}^n (x_k - \\bar x_n)=\\mathrm{var}(x)\n",
- "\\label{eq:var_estimate_i_think} \\tag{15}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Now we can calculate an estimate of the error\n",
- "$\\mathrm{err}_X^{\\phantom X}$ of the sample mean $\\bar x_n$:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{err}_X^2\n",
- "=\\frac{1}{n^2}\\sum_{ij} \\mathrm{cov}(X_i, X_j) \\nonumber\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\approx\\frac{1}{n^2}\\sum_{ij}\\frac{1}{n}\\mathrm{cov}(x) =\\frac{1}{n^2}n^2\\frac{1}{n}\\mathrm{cov}(x)\\nonumber\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\frac{1}{n}\\mathrm{cov}(x)\n",
- "\\label{eq:error_estimate} \\tag{16}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which is nothing but the sample covariance divided by the number of\n",
- "measurements in the sample.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, uncorrelated results\n",
- "\n",
- "In the special case that the measurements of the sample are\n",
- "uncorrelated (equivalently the stochastic variables $X_i$ are\n",
- "uncorrelated) we have that the off-diagonal elements of the covariance\n",
- "are zero. This gives the following estimate of the sample error:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{err}_X^2=\\frac{1}{n^2}\\sum_{ij} \\mathrm{cov}(X_i, X_j) =\n",
- "\\frac{1}{n^2} \\sum_i \\mathrm{var}(X_i),\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "resulting in"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\mathrm{err}_X^2\\approx \\frac{1}{n^2} \\sum_i \\mathrm{var}(x)= \\frac{1}{n}\\mathrm{var}(x)\n",
- "\\label{eq:error_estimate_uncorrel} \\tag{17}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "where in the second step we have used Eq. ([15](#eq:var_estimate_i_think)).\n",
- "The error of the sample is then just its standard deviation divided by\n",
- "the square root of the number of measurements the sample contains.\n",
- "This is a very useful formula which is easy to compute. It acts as a\n",
- "first approximation to the error, but in numerical experiments, we\n",
- "cannot overlook the always present correlations.\n",
- "\n",
- "\n",
- "\n",
- "## Statistics, computations\n",
- "For computational purposes one usually splits up the estimate of\n",
- "$\\mathrm{err}_X^2$, given by Eq. ([16](#eq:error_estimate)), into two\n",
- "parts"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "\\mathrm{err}_X^2 = \\frac{1}{n}\\mathrm{var}(x) + \\frac{1}{n}(\\mathrm{cov}(x)-\\mathrm{var}(x)),\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "which equals"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\frac{1}{n^2}\\sum_{k=1}^n (x_k - \\bar x_n)^2 +\\frac{2}{n^2}\\sum_{k\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation} \n",
- "=\\frac{\\tau}{n}\\cdot\\mathrm{var}(x)\n",
- "\\label{eq:error_estimate_corr_time} \\tag{19}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "and we see that $\\mathrm{err}_X$ can be expressed in terms the\n",
- "uncorrelated sample variance times a correction factor $\\tau$ which\n",
- "accounts for the correlation between measurements. We call this\n",
- "correction factor the *autocorrelation time*:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "$$\n",
- "\\begin{equation}\n",
- "\\tau = 1+2\\sum_{d=1}^{n-1}\\kappa_d\n",
- "\\label{eq:autocorrelation_time} \\tag{20}\n",
- "\\end{equation}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Statistics, effective number of correlations\n",
- "For a correlation free experiment, $\\tau$\n",
- "equals 1. From the point of view of\n",
- "eq. ([19](#eq:error_estimate_corr_time)) we can interpret a sequential\n",
- "correlation as an effective reduction of the number of measurements by\n",
- "a factor $\\tau$. The effective number of measurements becomes:"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "$$\n",
- "n_\\mathrm{eff} = \\frac{n}{\\tau}\n",
- "$$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "To neglect the autocorrelation time $\\tau$ will always cause our\n",
- "simple uncorrelated estimate of $\\mathrm{err}_X^2\\approx \\mathrm{var}(x)/n$ to\n",
- "be less than the true sample error. The estimate of the error will be\n",
- "too *good*. On the other hand, the calculation of the full\n",
- "autocorrelation time poses an efficiency problem if the set of\n",
- "measurements is very large.\n",
- "\n",
- "\n",
- "\n",
- "\n",
"\n",
"\n",
"## Linking the regression analysis with a statistical interpretation\n",
"\n",
- "Finally, we are going to discuss several statistical properties which can be obtained in terms of analytical expressions. \n",
+ "\n",
"The\n",
"advantage of doing linear regression is that we actually end up with\n",
"analytical expressions for several statistical quantities. \n",
@@ -5204,8 +3866,10 @@
},
{
"cell_type": "code",
- "execution_count": 30,
- "metadata": {},
+ "execution_count": 28,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from numpy import *\n",
@@ -5345,8 +4009,10 @@
},
{
"cell_type": "code",
- "execution_count": 31,
- "metadata": {},
+ "execution_count": 29,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"from numpy import *\n",
@@ -5491,8 +4157,10 @@
},
{
"cell_type": "code",
- "execution_count": 32,
- "metadata": {},
+ "execution_count": 30,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -5720,8 +4388,10 @@
},
{
"cell_type": "code",
- "execution_count": 33,
- "metadata": {},
+ "execution_count": 31,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
@@ -5789,8 +4459,10 @@
},
{
"cell_type": "code",
- "execution_count": 34,
- "metadata": {},
+ "execution_count": 32,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
@@ -5888,8 +4560,10 @@
},
{
"cell_type": "code",
- "execution_count": 35,
- "metadata": {},
+ "execution_count": 33,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"\"\"\"\n",
@@ -5974,8 +4648,10 @@
},
{
"cell_type": "code",
- "execution_count": 36,
- "metadata": {},
+ "execution_count": 34,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Common imports\n",
@@ -6068,8 +4744,10 @@
},
{
"cell_type": "code",
- "execution_count": 37,
- "metadata": {},
+ "execution_count": 35,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"# Common imports\n",
@@ -6150,8 +4828,10 @@
},
{
"cell_type": "code",
- "execution_count": 38,
- "metadata": {},
+ "execution_count": 36,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -6206,12 +4886,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" H = -J \\sum_{k}^L s_k s_{k + 1},\n",
- "\\label{_auto11} \\tag{21}\n",
+ "\\label{_auto2} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
@@ -6231,8 +4911,10 @@
},
{
"cell_type": "code",
- "execution_count": 39,
- "metadata": {},
+ "execution_count": 37,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -6277,12 +4959,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" H = - \\sum_j^L \\sum_k^L s_j s_k J_{jk}.\n",
- "\\label{_auto12} \\tag{22}\n",
+ "\\label{_auto3} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
@@ -6301,12 +4983,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\boldsymbol{H} = \\boldsymbol{X} J,\n",
- "\\label{_auto13} \\tag{23}\n",
+ "\\label{_auto4} \\tag{4}\n",
"\\end{equation}\n",
"$$"
]
@@ -6325,12 +5007,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\boldsymbol{y} = \\boldsymbol{X}\\boldsymbol{\\beta} + \\boldsymbol{\\epsilon},\n",
- "\\label{_auto14} \\tag{24}\n",
+ "\\label{_auto5} \\tag{5}\n",
"\\end{equation}\n",
"$$"
]
@@ -6344,8 +5026,10 @@
},
{
"cell_type": "code",
- "execution_count": 40,
- "metadata": {},
+ "execution_count": 38,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X = np.zeros((n, L ** 2))\n",
@@ -6369,12 +5053,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" C(\\boldsymbol{X}, \\boldsymbol{\\beta})= \\frac{1}{n}\\left\\{(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})^T(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})\\right\\}.\n",
- "\\label{_auto15} \\tag{25}\n",
+ "\\label{_auto6} \\tag{6}\n",
"\\end{equation}\n",
"$$"
]
@@ -6408,8 +5092,10 @@
},
{
"cell_type": "code",
- "execution_count": 41,
- "metadata": {},
+ "execution_count": 39,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X_train_own = np.concatenate(\n",
@@ -6424,8 +5110,10 @@
},
{
"cell_type": "code",
- "execution_count": 42,
- "metadata": {},
+ "execution_count": 40,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
@@ -6485,12 +5173,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\boldsymbol{\\beta} = \\boldsymbol{V}\\boldsymbol{\\Sigma}^{+} \\boldsymbol{U}^T \\boldsymbol{y}.\n",
- "\\label{_auto16} \\tag{26}\n",
+ "\\label{_auto7} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
@@ -6509,8 +5197,10 @@
},
{
"cell_type": "code",
- "execution_count": 43,
- "metadata": {},
+ "execution_count": 41,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
@@ -6520,8 +5210,10 @@
},
{
"cell_type": "code",
- "execution_count": 44,
- "metadata": {},
+ "execution_count": 42,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"beta = ols_svd(X_train_own,y_train)"
@@ -6536,8 +5228,10 @@
},
{
"cell_type": "code",
- "execution_count": 45,
- "metadata": {},
+ "execution_count": 43,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"J = beta[1:].reshape(L, L)"
@@ -6552,8 +5246,10 @@
},
{
"cell_type": "code",
- "execution_count": 46,
- "metadata": {},
+ "execution_count": 44,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -6597,12 +5293,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" H = -J \\sum_{k}^L s_k s_{k + 1},\n",
- "\\label{_auto17} \\tag{27}\n",
+ "\\label{_auto8} \\tag{8}\n",
"\\end{equation}\n",
"$$"
]
@@ -6618,8 +5314,10 @@
},
{
"cell_type": "code",
- "execution_count": 47,
- "metadata": {},
+ "execution_count": 45,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"import numpy as np\n",
@@ -6657,12 +5355,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" H = - \\sum_j^L \\sum_k^L s_j s_k J_{jk}.\n",
- "\\label{_auto18} \\tag{28}\n",
+ "\\label{_auto9} \\tag{9}\n",
"\\end{equation}\n",
"$$"
]
@@ -6681,12 +5379,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" H = X J,\n",
- "\\label{_auto19} \\tag{29}\n",
+ "\\label{_auto10} \\tag{10}\n",
"\\end{equation}\n",
"$$"
]
@@ -6705,12 +5403,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\boldsymbol{y} = \\boldsymbol{X}\\boldsymbol{\\beta} + \\boldsymbol{\\epsilon}.\n",
- "\\label{_auto20} \\tag{30}\n",
+ "\\label{_auto11} \\tag{11}\n",
"\\end{equation}\n",
"$$"
]
@@ -6724,8 +5422,10 @@
},
{
"cell_type": "code",
- "execution_count": 48,
- "metadata": {},
+ "execution_count": 46,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"X = np.zeros((n, L ** 2))\n",
@@ -6754,8 +5454,10 @@
},
{
"cell_type": "code",
- "execution_count": 49,
- "metadata": {},
+ "execution_count": 47,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"clf = skl.LinearRegression().fit(X_train, y_train)"
@@ -6770,8 +5472,10 @@
},
{
"cell_type": "code",
- "execution_count": 50,
- "metadata": {},
+ "execution_count": 48,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"J_sk = clf.coef_.reshape(L, L)"
@@ -6786,8 +5490,10 @@
},
{
"cell_type": "code",
- "execution_count": 51,
- "metadata": {},
+ "execution_count": 49,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -6820,8 +5526,8 @@
"metadata": {},
"source": [
"1\n",
- "9\n",
- "2\n",
+ "3\n",
+ "6\n",
" \n",
"<\n",
"<\n",
@@ -6842,8 +5548,10 @@
},
{
"cell_type": "code",
- "execution_count": 52,
- "metadata": {},
+ "execution_count": 50,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"_lambda = 0.1\n",
@@ -6874,12 +5582,12 @@
"metadata": {},
"source": [
"\n",
- "\n",
+ "\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" C(\\boldsymbol{X}, \\boldsymbol{\\beta}; \\lambda) = (\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y})^T(\\boldsymbol{X}\\boldsymbol{\\beta} - \\boldsymbol{y}) + \\lambda \\sqrt{\\boldsymbol{\\beta}^T\\boldsymbol{\\beta}}.\n",
- "\\label{_auto22} \\tag{32}\n",
+ "\\label{_auto13} \\tag{13}\n",
"\\end{equation}\n",
"$$"
]
@@ -6893,8 +5601,10 @@
},
{
"cell_type": "code",
- "execution_count": 53,
- "metadata": {},
+ "execution_count": 51,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"clf_lasso = skl.Lasso(alpha=_lambda).fit(X_train, y_train)\n",
@@ -6927,8 +5637,10 @@
},
{
"cell_type": "code",
- "execution_count": 54,
- "metadata": {},
+ "execution_count": 52,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"lambdas = np.logspace(-4, 5, 10)\n",
@@ -6990,8 +5702,10 @@
},
{
"cell_type": "code",
- "execution_count": 55,
- "metadata": {},
+ "execution_count": 53,
+ "metadata": {
+ "collapsed": false
+ },
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20, 14))\n",
@@ -7036,25 +5750,7 @@
]
}
],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "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.8.3"
- }
- },
+ "metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}