-
@@ -193,7 +215,7 @@ MathJax.Hub.Config({
+
-Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable. -
-Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as -$$ -\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right), -$$ + +
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
-where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
-
-
-
@@ -193,6 +222,11 @@ where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as
- +
-In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring -$$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0, -$$ -which results in -$$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0, -$$ - -or in a matrix-vector form as -$$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right). -$$ - -where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \). -
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[2]])
+ypredict = linreg.predict(xnew)
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
@@ -197,6 +219,12 @@ where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix
+
-We can rewrite -$$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right), -$$ - -as -$$ -\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta}, -$$ - -and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution -$$ -\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}. -$$ -
@@ -194,6 +202,11 @@ $$
-If we then introduce the matrix +Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable. + +
+Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as $$ -\hat{H} = \hat{A}^T\hat{A}, +\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right), $$ -we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{H} \) are \( h_{ij} \)) -$$ -\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik} -$$ +where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements. -We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as -$$ -\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2, -$$ - -resulting in -$$ -\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}! -$$ +
@@ -198,6 +212,11 @@ $$
-The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write + +
+In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring $$ -y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i. +\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0, $$ -By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by +which results in $$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0, +\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0, $$ -and +or in a matrix-vector form as $$ -\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0. +\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right). $$ + +where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \).
-We define then +We can rewrite $$ -\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2}, +\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right), $$ - +as $$ -\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2}, +\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta}, $$ +and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution $$ -\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right), +\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}. $$ - -$$ -\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2}, -$$ - -$$ -\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2}, -$$ - -and show that -$$ -\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}, -$$ - -$$ -\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}. -$$ - -
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed below. @@ -210,6 +213,11 @@ The LSM suffers often from both being underdetermined and overdetermined in the
-How can we use the singular value decomposition to find the parameters \( \beta_j \)? More details will come. We first note that a general \( m\times n \) matrix \( \hat{A} \) can be written in terms of a diagonal matrix \( \hat{\Sigma} \) of dimensionality \( n\times n \) and two orthognal matrices \( \hat{U} \) and \( \hat{V} \), where the first has dimensionality \( m \times n \) and the last dimensionality \( n\times n \). We have then + +
+If we then introduce the matrix $$ -\hat{A} = \hat{U}\hat{\Sigma}\hat{V} +\hat{H} = \hat{A}^T\hat{A}, +$$ + +we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{H} \) are \( h_{ij} \)) +$$ +\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik} +$$ + +We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as +$$ +\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2, +$$ + +resulting in +$$ +\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}! $$
-
diff --git a/doc/pub/Regression/html/Regression-bs.html b/doc/pub/Regression/html/Regression-bs.html index cc3f0021f..4a67b29fe 100644 --- a/doc/pub/Regression/html/Regression-bs.html +++ b/doc/pub/Regression/html/Regression-bs.html @@ -73,13 +73,30 @@ Automatically generated HTML file from DocOnce source 2, None, '___sec10'), - ('The $\\chi^2$ function', 2, None, '___sec11'), - ('The $\\chi^2$ function', 2, None, '___sec12'), - ('The $\\chi^2$ function', 2, None, '___sec13'), + ('Simple regression model', 2, None, '___sec11'), + ('Simple regression model, now using _scikit-learn_', + 2, + None, + '___sec12'), + ('Correlations and the quality of our results', + 2, + None, + '___sec13'), ('The $\\chi^2$ function', 2, None, '___sec14'), ('The $\\chi^2$ function', 2, None, '___sec15'), ('The $\\chi^2$ function', 2, None, '___sec16'), - ('The singular value decompostion', 2, None, '___sec17')]} + ('The $\\chi^2$ function', 2, None, '___sec17'), + ('The $\\chi^2$ function', 2, None, '___sec18'), + ('The $\\chi^2$ function', 2, None, '___sec19'), + ('Simple regression model with gradient descent', + 2, + None, + '___sec20'), + ('Simple regression model with stochastic gradient descent', + 2, + None, + '___sec21'), + ('The singular value decompostion', 2, None, '___sec22')]} end of tocinfo --> @@ -128,13 +145,18 @@ MathJax.Hub.Config({
-
@@ -193,7 +215,7 @@ MathJax.Hub.Config({
-
+ + +
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Linear Regression')
+plt.show()
++ + +
# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[2]])
+ypredict = linreg.predict(xnew)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
+$$
+\mathrm{Corr}(X,Y) = \frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2} }
+$$
+
+
+
+Another quantity is the autocorrelation function we discussed in our chapter on statistical analysis. +
@@ -515,7 +592,7 @@ where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as
@@ -546,7 +623,7 @@ where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix
@@ -575,7 +652,7 @@ $$
@@ -611,7 +688,7 @@ $$
@@ -640,7 +717,7 @@ $$
@@ -695,7 +772,78 @@ The LSM suffers often from both being underdetermined and overdetermined in the
+
+
+
+
+
+
diff --git a/doc/pub/Regression/html/Regression-solarized.html b/doc/pub/Regression/html/Regression-solarized.html
index a0ba5efaa..154f6eca6 100644
--- a/doc/pub/Regression/html/Regression-solarized.html
+++ b/doc/pub/Regression/html/Regression-solarized.html
@@ -93,13 +93,30 @@ div { text-align: justify; text-justify: inter-word; }
2,
None,
'___sec10'),
- ('The $\\chi^2$ function', 2, None, '___sec11'),
- ('The $\\chi^2$ function', 2, None, '___sec12'),
- ('The $\\chi^2$ function', 2, None, '___sec13'),
+ ('Simple regression model', 2, None, '___sec11'),
+ ('Simple regression model, now using _scikit-learn_',
+ 2,
+ None,
+ '___sec12'),
+ ('Correlations and the quality of our results',
+ 2,
+ None,
+ '___sec13'),
('The $\\chi^2$ function', 2, None, '___sec14'),
('The $\\chi^2$ function', 2, None, '___sec15'),
('The $\\chi^2$ function', 2, None, '___sec16'),
- ('The singular value decompostion', 2, None, '___sec17')]}
+ ('The $\\chi^2$ function', 2, None, '___sec17'),
+ ('The $\\chi^2$ function', 2, None, '___sec18'),
+ ('The $\\chi^2$ function', 2, None, '___sec19'),
+ ('Simple regression model with gradient descent',
+ 2,
+ None,
+ '___sec20'),
+ ('Simple regression model with stochastic gradient descent',
+ 2,
+ None,
+ '___sec21'),
+ ('The singular value decompostion', 2, None, '___sec22')]}
end of tocinfo -->
-
+
+
+
+
+
+
+
+
+Another quantity is the autocorrelation function we discussed in our chapter on statistical analysis.
+
+
+
@@ -464,7 +554,7 @@ where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as
@@ -492,7 +582,7 @@ where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix
@@ -518,7 +608,7 @@ $$
@@ -549,7 +639,7 @@ $$
@@ -573,7 +663,7 @@ $$
@@ -618,7 +708,76 @@ The LSM suffers often from both being underdetermined and overdetermined in the
+
+
+
+
+
+
+
+
@@ -635,7 +794,7 @@ $$
-
+
+
+
+
+
+
+
+
+Another quantity is the autocorrelation function we discussed in our chapter on statistical analysis.
+
+
+
@@ -469,7 +559,7 @@ where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as
@@ -497,7 +587,7 @@ where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix
@@ -523,7 +613,7 @@ $$
@@ -554,7 +644,7 @@ $$
@@ -578,7 +668,7 @@ $$
@@ -623,7 +713,76 @@ The LSM suffers often from both being underdetermined and overdetermined in the
+
+
+
+
+
+
+
+
@@ -640,7 +799,7 @@ $$
The \( \chi^2 \) function
+The \( \chi^2 \) function
The \( \chi^2 \) function
+The \( \chi^2 \) function
The \( \chi^2 \) function
+The \( \chi^2 \) function
The \( \chi^2 \) function
+The \( \chi^2 \) function
The \( \chi^2 \) function
+The \( \chi^2 \) function
The singular value decompostion
+Simple regression model with gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+theta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 1000
+m = 100
+
+for iter in range(Niterations):
+ gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
+ theta -= eta*gradients
+
+print(theta)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+ypredict2 = xbnew.dot(theta_linreg)
+plt.plot(xnew, ypredict, "r-")
+plt.plot(xnew, ypredict2, "b-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
Simple regression model with stochastic gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import SGDRegressor
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
+sgdreg.fit(x,y.ravel())
+print(sgdreg.intercept_, sgdreg.coef_)
+
The singular value decompostion
Nov 26, 2017
May 14, 2018
@@ -441,7 +458,80 @@ meaning that the solution for \( \hat{\beta} \) is the one which minimizes the r
-The \( \chi^2 \) function
+Simple regression model
+We are now ready to write our first program which aims at solving the above linear regression equations. We start with data we have produced ourselves, in this case normally distributed random numbers along the \( x \)-axis. These numbers define then the value of a function \( y(x)=4+3x+N(0,1) \). Thereafter we order the \( x \) values and employ our linear regression algorithm to set up the best fit.
+
+# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Linear Regression')
+plt.show()
+
+
+Simple regression model, now using scikit-learn
+We can repeat the above algorithm using scikit-learn as follows.
+# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[2]])
+ypredict = linreg.predict(xnew)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
+
+Correlations and the quality of our results
+In order to test the quality of our fit, there are several measures which can be implemented. One is the so-called
+correlation function defined as
+$$
+\mathrm{Corr}(X,Y) = \frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2} }
+$$
+
+
+
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The singular value decompostion
+Simple regression model with gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+theta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 1000
+m = 100
+
+for iter in range(Niterations):
+ gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
+ theta -= eta*gradients
+
+print(theta)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+ypredict2 = xbnew.dot(theta_linreg)
+plt.plot(xnew, ypredict, "r-")
+plt.plot(xnew, ypredict2, "b-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
+
+Simple regression model with stochastic gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import SGDRegressor
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
+sgdreg.fit(x,y.ravel())
+print(sgdreg.intercept_, sgdreg.coef_)
+
+
+The singular value decompostion
Nov 26, 2017
May 14, 2018
@@ -446,7 +463,80 @@ meaning that the solution for \( \hat{\beta} \) is the one which minimizes the r
-The \( \chi^2 \) function
+Simple regression model
+We are now ready to write our first program which aims at solving the above linear regression equations. We start with data we have produced ourselves, in this case normally distributed random numbers along the \( x \)-axis. These numbers define then the value of a function \( y(x)=4+3x+N(0,1) \). Thereafter we order the \( x \) values and employ our linear regression algorithm to set up the best fit.
+
+# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Linear Regression')
+plt.show()
+
+
+Simple regression model, now using scikit-learn
+We can repeat the above algorithm using scikit-learn as follows.
+# Importing various packages
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import LinearRegression
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+linreg = LinearRegression()
+linreg.fit(x,y)
+xnew = np.array([[0],[2]])
+ypredict = linreg.predict(xnew)
+
+plt.plot(xnew, ypredict, "r-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
+
+Correlations and the quality of our results
+In order to test the quality of our fit, there are several measures which can be implemented. One is the so-called
+correlation function defined as
+$$
+\mathrm{Corr}(X,Y) = \frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2} }
+$$
+
+
+
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The \( \chi^2 \) function
+The \( \chi^2 \) function
-The singular value decompostion
+Simple regression model with gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+theta = np.random.randn(2,1)
+
+eta = 0.1
+Niterations = 1000
+m = 100
+
+for iter in range(Niterations):
+ gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
+ theta -= eta*gradients
+
+print(theta)
+xnew = np.array([[0],[2]])
+xbnew = np.c_[np.ones((2,1)), xnew]
+ypredict = xbnew.dot(theta)
+ypredict2 = xbnew.dot(theta_linreg)
+plt.plot(xnew, ypredict, "r-")
+plt.plot(xnew, ypredict2, "b-")
+plt.plot(x, y ,'ro')
+plt.axis([0,2.0,0, 15.0])
+plt.xlabel(r'$x$')
+plt.ylabel(r'$y$')
+plt.title(r'Random numbers ')
+plt.show()
+
+
+Simple regression model with stochastic gradient descent
+Add info about the equations, play around with different learning rates
+# Importing various packages
+from math import exp, sqrt
+from random import random, seed
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn.linear_model import SGDRegressor
+
+x = 2*np.random.rand(100,1)
+y = 4+3*x+np.random.randn(100,1)
+
+xb = np.c_[np.ones((100,1)), x]
+theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
+print(theta_linreg)
+sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
+sgdreg.fit(x,y.ravel())
+print(sgdreg.intercept_, sgdreg.coef_)
+
+
+The singular value decompostion