<jupyter_start><jupyter_text>Linear regressionGenerate a random data with underlying trend: $y = 4 + 3x$.By convention, the y-intercept is taken to be the first parameter.<jupyter_code>import numpy as np
X= 2* np.random.rand(100, 1)
y = 4 + 3*X + np.random.randn(100, 1)
import matplotlib.pyplot as plt
plt.scatter(X,y)<jupyter_output><empty_output><jupyter_text>Our model is:  $y = c + mx = \theta^T X$ where $\theta = (c, m)$ and $X =  (1, x)$.The given data is$$ y_i = \theta^T X_i$$for $i = 1, ..., M$.We find the the optimum parameter $\hat{\theta}$ for which $$ MSE(\theta) = \frac1m \sum_{i=1}^{m} \left( \theta^T . X_i - y_i \right)^2$$is minimized.Observe that$$\begin{bmatrix}\theta^T X_1 \\\vdots \\\theta^T X_M \\ \end{bmatrix} -\begin{bmatrix}y_1\\\vdots \\y_M \end{bmatrix}  = \bar{X}\theta - y. $$ where $$ \bar{X} =  \begin{bmatrix} X_1\\ \vdots\\ X_M \end{bmatrix}=  \begin{bmatrix} 1 &x_1\\ \vdots & \vdots\\1 & x_m \end{bmatrix}. $$So,$$\text{MSE}(\theta) = \frac1m (\bar{X}\theta - y)^T(\bar{X}\theta - y)$$Simplifying,$$ M. \text{MSE}(\theta) = (\bar{X}\theta - y)^T(\bar{X}\theta - y) = (\bar{X}\theta)^T (\bar{X} \theta) - 2Y^T \bar{X}^T + Y^T Y $$Using the fact that $$ \frac{d}{dX} X^TX = 2X$$and $$ \frac{d}{dX} X^TB = B$$we set  the derivative of MSE to zero to get$$ 2 \bar{X}^T X \theta = 2 \bar{X}^T Y$$which yields the normal equation$$\hat{\theta} = \left(\bar{X}^T \bar{X} \right)^{-1} \bar{X}^T y$$<jupyter_code># Adding bias column of ones, the last column of \bar{X}
X_b = np.c_[np.ones((100,1)), X] # Note the double (( )) !!
# using numpy function linalg.inv
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
theta_best<jupyter_output><empty_output><jupyter_text>Our expection was 4 and 3. Now, we compute the prediction of the model, which is a line.<jupyter_code>y_predict = np.c_[np.ones((100,1)), X].dot(theta_best)
plt.scatter(X,y)
plt.plot(X, y_predict,'-r')<jupyter_output><empty_output><jupyter_text>In `scikit-learn`, we can do the same using `LinearRegression`. We feed the x-values without the bias column and `lin_reg.intercept_` outputs the intercept.While predicting, we need to take the dot product with x-value with a bias column (values of 1).<jupyter_code>from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X,y)<jupyter_output><empty_output><jupyter_text>**Note:** You do *not* have to manually add the bias column to use `LinearRegression()`. The fitted parameters can be viewed using the variable with underscore in the end.<jupyter_code>lin_reg.intercept_, lin_reg.coef_<jupyter_output><empty_output><jupyter_text>Gradient descentThe **learning rate**  $\eta$ is the step size by which the gradient descent method jumps to the next iteration. The gradient descent algorithm is$$ \theta^{n+1} = \theta^n - \eta \nabla_{\theta} \text{MSE}(\theta^n)$$where$$ \nabla_\theta \text{MSE}(\theta) = \frac2M  \bar{X}^T(\bar{X}\theta -Y)$$The inital guess is chosen randomly. However, depending on the guess, might take a long time to reach the minimizer.The **tolerance** is the lower bound of the gradient of MSE.Here, `X_b`=$\bar{X}$.<jupyter_code>eta = 0.1
n_iter = 1000
m = 100
theta = np.random.randn(2,1)

for iteration in range(n_iter):
    gradient = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta  - eta*gradient

theta<jupyter_output><empty_output><jupyter_text>We can do a grid search to find the best learning rate that works. (Not sure how that will give us a different value of the minimizer, but I'm guessing it will be better for MSE's with multiple local minimizers.)  For a convex cost function, (smooth enough), for a fixed learning rate, the convergenge rate is $O(\frac{1}{\\ iterations})$. Stochastic gradient descentHere, the elements of the gradient (a vectot) is computed using $$ \frac{\partial}{\partial\theta_j} \text{MSE}(\theta) = \frac2M \sum_{i=1}^m (\theta^T X_i - y_i) X_i(j)$$for each $j$. But for each $j$, a different sample is chosen. (?)Choose the gradient of MSE from randomly chosen samples. However, we need to reduce the learning rate gradually to make sure randomly chosen gradients do not make the guess run out of the theorerical minimizer too much.This method does not converge to the minimizer as step size goes to infinity, but it *does* converge in an average.On the other hand, exploiting the randomness, it is very fast. **Note:** `.ravel()` is same as `.reshape(-1)` which is used to convert a matrix (with ``[[]]``) to an 1-dimensional array (with `[]`) which is required for `SGDRegressor()`.<jupyter_code>from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(tol=0.2, penalty=None, eta0=0.1)
sgd_reg.fit(X, y.ravel())
sgd_reg.intercept_, sgd_reg.coef_<jupyter_output><empty_output><jupyter_text>Mini Batch gradient descentchooses a subset of the whole dataset to compute the gradient. This dataset is chosen at random so it has the best of both Stochastic and Batch gradient descent (where the gradient is computed from the whole dataset), fast and converges to the minimizer.  Polynomial regressionLet the model be$$ y = 0.5x^2 + x+ 2$$<jupyter_code>m = 100
X = 6 * np.random.rand(m, 1) - 3
y = 0.5 * X**2 + X + 2 + np.random.randn(m,1)
plt.scatter(X,y)<jupyter_output><empty_output><jupyter_text>Here, we create columns including all higher powers of $x$ upto the `degree` provided, with or without the bias term (i.e. the constant coefficent).For example, after fitting the `PolynomialFeatures()`, the first column (with `include_bias=False`) is $x$, second column is $x^2$ etc.<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=2, include_bias = False)
X_poly = poly_features.fit_transform(X)
X_poly[0]  # contains x and x^2<jupyter_output><empty_output><jupyter_text>A polynomial regression is essentially a linear regression with more multiple features where the other features are related to one feature in a polyminal way.<jupyter_code>lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
lin_reg.intercept_, lin_reg.coef_<jupyter_output><empty_output><jupyter_text>This implies that the model is $y = 1.9 + 0.92 x + 0.52 x^2$. If $x=(a,b)$ is two-dimensionals, `PolynomialFeatures` with `degree=2` will produce $a, b, a2, b^2, ab$ columns. For degree $d$ with $n$ features, total polynomial features will be$$ \frac{(n+d)!}{d! n!}$$ features. **Note:** Just like the linear case (i.e. when your model or hypothesis is linear), you do *not* have to add the bias column to your data manually to use the `LinearRegression()`. We plot the data and the fitting polynomial.<jupyter_code># equally space 100 numbers between the min and max of the data
X_arr = np.linspace(min(X), max(X), num=100)
# collecting the coefficients of polynomial
reg_params = np.c_[lin_reg.intercept_, lin_reg.coef_]
# computing the y-values of the polynomial
y_curve = np.c_[np.ones((100,1)), X_arr, np.power(X_arr,2)].dot(np.transpose(reg_params))
plt.plot(X_arr, y_curve, 'r')
plt.scatter(X, y)<jupyter_output><empty_output><jupyter_text>Alternate way to plot the same curve: by predicting an equally-spaced 100 `linspace()` values using the fitted model.<jupyter_code>X_arr = np.linspace(min(X), max(X), num=100)
X_poly_arr = poly_features.fit_transform(X_arr)
plt.plot(X_arr, lin_reg.predict(X_poly_arr), 'r')
plt.scatter(X, y)<jupyter_output><empty_output><jupyter_text>Ploynomial regression using a Stochastic gradient descentWe can do the polynomial regression using a Stochastic Gradient Descent since a polynomial regression is essentially a *linear regression* with multiple columns (where the columns were constructed by computing powers of the main column).<jupyter_code>sgd_poly_reg = SGDRegressor(tol=0.2, penalty=None, eta0=0.1)
sgd_poly_reg.fit(X_poly,y.reshape(-1)) # SGDRegressor requires y to be a flattened array, not a vector
plt.plot(X_arr, sgd_poly_reg.predict(X_poly_arr))
plt.scatter(X, y)<jupyter_output><empty_output><jupyter_text>Note: every time we run it, the curve it going to be slightly different since we are using *stochastic* method.  How to determine that the model is under or over fitting?* Look at the cross-validation scores. If the model fits a small amount of data *very* well (ie. interpolating), it would not work well with randomly selected data (overfitting).1. **Learning curves:** Plot the graph of error of the training and test data vs the sample size of the set.<jupyter_code>from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

def plot_learning_curves(model, X, y):
    X_train, X_validation, y_train, y_validation = train_test_split(X, y, test_size=0.2)
    train_errors, validation_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        training_prediction = model.predict(X_train[:m])
        validation_prediction = model.predict(X_validation)
        train_errors.append(mean_squared_error(training_prediction, y_train[:m]))
        validation_errors.append(mean_squared_error(validation_prediction, y_validation))
    plt.plot(np.sqrt(train_errors), "r-+", label="train" )
    plt.plot(np.sqrt(validation_errors), "g-", label="validation" )
lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)<jupyter_output><empty_output><jupyter_text>Since the error from the training data is less than the error from the validation data (there is a gap for all size of training dataset), we are underfitting.When the training dataset is small, we are overfitting, hence the error on the training dataset is very small. At this stage, the model hasn't seen the validation set, so the error on that is very high. As the dataset size increases, the error catches up.<jupyter_code>from sklearn.pipeline import Pipeline
polynomial_regression = Pipeline([
    ("poly_features", PolynomialFeatures(degree=10, include_bias=False)),
    ("sgd_reg", LinearRegression()),

])
plot_learning_curves(polynomial_regression, X, y)<jupyter_output><empty_output><jupyter_text>Plotting the corresponding models<jupyter_code>X_arr = np.linspace(min(X), max(X), num=100)
plt.plot(X_arr, lin_reg.predict(X_arr),'r')
plt.plot(X_arr, polynomial_regression.predict(X_arr),'b' )
plt.scatter(X, y)<jupyter_output><empty_output><jupyter_text>Ridge Regression or Tikhonov RegularizationThe cost function is$$ J(\theta) = \text{MSE}(\theta) + \frac\alpha2 \lvert \lvert \theta \rvert \rvert_{l^2}$$ Or,$$ J(\theta) = (\bar{X}\theta - y)^T ( \bar{X} \theta - y) + \frac\alpha2 \theta^T \theta$$ For $\alpha=0$, the method reduces to Linear Regression. As $\alpha$ increases, the fitting curve becomes flatter and eventually approaches a horizontal straight line.**Note:** Requires scaling before we apply Ridge regression.<jupyter_code>from sklearn.linear_model import Ridge
ridge_reg = Ridge(alpha=1)
ridge_reg.fit(X, y)
ridge_reg.predict([[1.5]])<jupyter_output><empty_output><jupyter_text>Lasso Regression or Selection Operator RegressionThe cost function is$$ J(\theta) = \text{MSE}(\theta) + \alpha \lvert \lvert \theta \rvert \rvert_{l^1}$$ Or,$$ J(\theta) = (\bar{X}\theta - y)^T ( \bar{X} \theta - y) + \alpha \sum_{i=1}^n \lvert \theta_i \rvert$$ The *subgradient* (since the gradient does not exist when $\theta_i = 0$ for some $i$) of the cost function is$$ g(\theta, J) = \nabla_\theta \text{MSE}(\theta) + \alpha\begin{bmatrix}sgn(\theta_1)\\\vdots \\sgn(\theta_n)\end{bmatrix}$$**Advantage:** The Lasso regression zeros out the smallest (least important) parameters. In other words, it does a **feature selection** and ouputs a **sparse model**.**Disadvatage:** Behaves badly when the number of features is bigger than the number of the datapoints, or when several features are strongly correlated.<jupyter_code>from sklearn.linear_model import Lasso
lasso_reg = Lasso(alpha=0.1)
lasso_reg.fit(X,y)
lasso_reg.predict([[1.5]])<jupyter_output><empty_output><jupyter_text>Elastic NetThe regularization involves a convex combination of the $l^1$ and $l^2$ norm of $\theta$. The cost fuction is$$ J(\theta) = \text{MSE}(\theta) + \alpha \left( r \lvert \lvert \theta \rvert \rvert_{l^2} + (1-r) \lvert \lvert \theta \rvert \rvert_{l^1} \right)$$for $r \in [0, 1]$. Usage of linear regression modelsUse Elastic Net or Lasso if there are reasons to suspect that only a few features are useful.**Ridge > Elastic Net > Lasso > Linear Regression**  Early Stopping as a regularization methodWe stop the learning process when the validation error reaches a local minimum. **Beautiful free lunch**However, since the error is not smooth, let the process run for a few more steps to make sure that the error does not go down again and roll back the parameters later.  Note sure what the next code does. How is this `epoch` being inside the loop?<jupyter_code>X_train, X_validation, y_train, y_validation = train_test_split(X, y, test_size=0.2)
X_train_poly = poly_features.fit_transform(X_train)

from sklearn.preprocessing import StandardScaler
std_sc = StandardScaler()
X_train_poly_scaled = std_sc.fit_transform(X_train_poly)

from sklearn.base import clone
sgd_reg = SGDRegressor(n_iter=1, warm_start=True, penalty=None,
                       learning_rate="constant", eta0=0.0005,
                      )
minimum_val_error = float("inf")
best_epoch = None
best_model = None

for epoch in range(1000):
    sgd_reg.fit(X_train_poly_scaled, y) #warm starts from last cloned point 
    y_val_predict = sgd_reg.predict(std_sc.fit_transform(poly_features.fit_transform(X_validation)))
    val_error = mean_squared_error(y_val_predict, y_validation)
    if val_error < minimum_val_error:
        minimum_val_error = val_error
        best_epoch = epoch
        best_model = clone(sgd_reg)<jupyter_output>/usr/local/lib/python3.6/dist-packages/sklearn/linear_model/stochastic_gradient.py:152: DeprecationWarning: n_iter parameter is deprecated in 0.19 and will be removed in 0.21. Use max_iter and tol instead.
  DeprecationWarning)
/usr/local/lib/python3.6/dist-packages/sklearn/utils/validation.py:761: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)<jupyter_text>Logistic Regression for classificationA classification problem can be thought as a regression problem where we fit a curve in the data and group the data points on each side of the curve as of a single class. A logistic curve.<jupyter_code>points = np.linspace(-10, 10)
plt.plot(points, np.power(1 + np.exp(-points), -1))<jupyter_output><empty_output><jupyter_text>Let the boundary of the 2-label classification (or the regression curve to fit the data) be $y_i = \theta^T X_i$ for each $i=1, \dots, m$.Define the probability $p_i$ that the point $(X_i, y_i)$ lies on the *positive* side of the regression curve (i.e. when $\theta^T X_i > 0$) by$$p_i = \left(1 + e^{-\theta^T X_i}\right)^{-1}.$$ When $p_i = \frac12$,  the point $(X_i, y_i)$ lies exactly on the boundary between the two classes (or, interpolated exactly by the model $y_i = \theta^T X_i$), since the range of the logistic funcion $(1 + e^{-t})^{-1}$ is $[0, 1]$ with the midpoint at $\frac12$.So, we can state our classification prediction algorithm as$$y_i = \begin{cases}0 \text{ if } p_i < 1/2\\1 \text{ if } p_i > 1/2,\end{cases}$$although we can just work with the actual probability $p_i$.The cost value for each point is defined by$$\begin{align*}c_i(\theta) & = \begin{cases}-\log(p_i) & \text{ if } y_i = 1\\ -\log(1 - p_i)  &\text{ if } y_i = 0 \end{cases}\\%& =%\begin{cases}%\log(1 + e^{\theta^T X_i})  & \text{ if }  (1 + e^{-\theta^T X_i})^{-1} > \frac12\\ %\log(1 + e^{-\theta^T X_i})  &\text{ if }  (1 + e^{-\theta^T X_i})^{-1} < \frac12.%\end{cases}\end{align*}$$This choice of cost value reflects that whenever $p_i \in [0, 1]$ is small (smaller than $\frac12$), i.e. when the datapoint lies on the *negative* side of the regression curve, the cost value corresponding to $y_i = 1$ (i.e. the point is on the *positive* side of the regression curve) is bigger and vice versa. The cost function (called *log loss*) is defined as the sum of the convex combination of the cost values of each datapoint$$J(\theta) = -\frac1m \sum_{i=1}^m \text{floor}\left( y_i \log(p_i) + (1 - y_i) \log( 1 - p_i) \right)$$The cost function is convex in $\theta$ hence (?) we can differentiate and find the gradient.$$\frac{\partial}{\partial \theta_j} J(\theta) = \frac1m \sum_{i=1}^m \left( \left( 1 + e^{-\theta^T X_i} \right)^{-1} - y_i \right) (X_i)_j$$Using this gradient, we can train the model with gradient descent or batch gradient descent method (with randomly chosed set of indices to compute the gradient) using$$\theta^{n+1} = \theta^n - \eta \frac{\partial}{\partial \theta_j} J(\theta^n)$$  In the notation of probabilityIn terms of probability, we have$$\mathbb{P}(y_i = 1 | X_i) = \left( 1 + e^{-\theta^T X_i}\right)^{-1}\\$$and$$\mathbb{P}(y_i = 0 | X_i) = 1 - \left( 1 + e^{-\theta^T X_i}\right)^{-1} = \frac{e^{-\theta^T X_i}}{1 + e^{-\theta^T X_i}}$$Generalizing, for any $y \in [0,1]$, we can write,$$\mathbb{P}(y_i = y | X_i) = \left( 1 + e^{-\theta^T X_i}\right)^{-y}\left( \frac{e^{-\theta^T X_i}}{1 + e^{-\theta^T X_i}}\right)^{(1-y)}$$The **likelihood function** (not sure about the definition)$$L = \prod_{i=1}^{m} \mathbb{P}(y_i = y_i| X_i)$$The **log-likelihood function** is $$\log L.$$ Using the Iris dataset.<jupyter_code>from sklearn import datasets
iris = datasets.load_iris()
list(iris.keys())
iris['feature_names']
# select petal width (cm)
X = iris['data'][:, 3:]
y = (iris['target'] == 2).astype(np.int)
y<jupyter_output><empty_output><jupyter_text>Train a logistic regression model to train type 2 flower by the petal width.<jupyter_code>from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression()
log_reg.fit(X,y)
# creating a uniformly spaced vector of petal widths
X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
# predict_proba instead of predict to output the proboblity instead of the class
y_proba = log_reg.predict_proba(X_new)
y_proba<jupyter_output><empty_output><jupyter_text>i-th row prints the probablities that the flower with petal width `X_new[i]` is a type-2 flower and not a type-2 flower.<jupyter_code>plt.plot(X_new, y_proba[:, 1], 'g', label='yes' )
plt.plot(X_new, y_proba[:, 0], 'b', label = 'no')<jupyter_output><empty_output><jupyter_text>Where the two graphs intersects (i.e. around 1.6), the decision boundary is set, i.e., higher than 1.6 petal width implies the model predicts  that the flower is indeed of type-2.  Now we train with more features<jupyter_code>log_reg2 = LogisticRegression()
log_reg2.fit(iris['data'][:, 0:3], y)
# predicted probability
print(log_reg2.predict_proba(iris['data'][50, 0:3].reshape(1,-1)))
print(iris['target'][50])  # actual label<jupyter_output>[[0.913103 0.086897]]
1<jupyter_text>Can do a regularized cost function. By defalult, the model is regularized by (`penalty=`) the `l2` norm of $\theta$ with the regularization constant `C` = $\frac1\alpha$. Higher `C` implies less regularized model.  Softmax Regressionfor classifying multiple labels To classify multiple ($K$) classes, a probability is computed from multiple linear regression models and the class with the maximum probability is returned.$$ S_k(X_i) = \theta_k^T X_i$$Probability that the given data sample $X_i$ belongs to the $k$th class is defined by $$ p_i^k = \frac{e^{s_k(X_i)}}{\sum_{j=1}^K e^{s_j(X_i)}}$$For predicting classes, we computer the class of the input by the class that produces the highest probability$$ \hat{y_i} = \text{argmax} \{ p_k^i \}$$The cost function called **cross-entropy** is given by$$J(\Theta) = - \frac1m \sum_{i=1}^{m} \sum_{k=1}^{K} y_i^k \log( p_i^k )$$where $$\Theta = (\theta_1, \cdots, \theta_K)$$and$$ y_i^k = \begin{cases}1 \text{ if } y_i \in \text{ Class} (k)\\0 \text{ otherwise}\end{cases}$$ The gradient of the cost function with respect to the parameter vector associated with the $k$-th class is given by$$\nabla_{\theta_k} J(\Theta) = \frac{1}{m} \sum_{i=1}^m \left( p_i^k - y_i^k\right) X_i$$  ImplementationThe softmax regression is `LogisgticRegression()` with parameter `multi_class="multinomial"` and `solver="lbftgs"`.By default, it adds $l^2$ regularization (i.e. $J(\Theta) + c ||\Theta||_2 $).We use *petal_length* and *petal_width*  (column 2 and 3) to train the softmax classifier to classify 3 types of flowers.<jupyter_code>X = iris['data'][:, (2,3)]
y= iris['target']

softmax_reg = LogisticRegression(multi_class='multinomial', solver='lbfgs', C=10)
softmax_reg.fit(X, y)<jupyter_output><empty_output><jupyter_text>If the *petal_length* and *petal_width* are 5 and 2 respectively, what kind of flower is it?<jupyter_code>softmax_reg.predict([[5,2]])<jupyter_output><empty_output><jupyter_text>What is the probability that it is of flower type 0, 1 and 2?<jupyter_code>softmax_reg.predict_proba([[5,2]])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1360822778782973, 'fraction_numerical': 0.015427102180744692, 'mean_word_length': 5.1161910308677925, 'pattern_counts': {'":': 0, '<': 110, '<?xml version=': 0, '>': 113, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '408052', 'license': 'MIT', 'n_tokens_mistral': 7407, 'n_tokens_neox': 7040, 'n_words': 2315, 'path': 'Models.ipynb', 'repo_name': 'debdeepbh/ml', 'size': '21002'}