<jupyter_start><jupyter_text>Regression is Prediction Linear RegressionYou all know the equation for a line, $y = mx+b$, or in the language of ML $f(x) = mx+b$. Two points are "produced" by $f(x)$:<jupyter_code># LAMBDA SCHOOL
# 
# MACHINE LEARNING
#
# MIT LICENSE

import numpy as np
import matplotlib.pyplot as plt

x = [3,8]
y = [4,7]

def plot_two(x, y):

    plt.figure(figsize=(3,2))
    plt.plot(x,y,'o')
    plt.xlim(0,10)
    plt.ylim(0,10)

plot_two(x, y)<jupyter_output><empty_output><jupyter_text>It is easy to solve $f(x) = mx+b$ for the above two points. We don't have to! There are many implementations of linear regression in the `python` universe.<jupyter_code>p = np.polyfit(x,y,1)<jupyter_output><empty_output><jupyter_text>So, $f(x) = \frac{3}{5}x + 2\frac{1}{5}$<jupyter_code># Really weird - the line drawn in this graph doesn't quite intersect the points
# how I expected. I tweaked the values to make it suit my Mark-One-Eyeball
# better. I suspect that the discrepancy is actually caused by the graph.
def f(x):
    return 0.595*x + 2.16
plot_two(x, y)
xs = np.linspace(0,10,11)
plt.plot(xs,f(xs),c='g');<jupyter_output><empty_output><jupyter_text>Co-LinearityThe above two points we used to generate $f(x)$ lie on the same line. Any pair of points lie on the same line - a line is defined by two points. Every other point on the line is co-linear with the first two. If we had three points all defined by $f(x)$, the result of `np.polyfit` would be the same.What about when the points are not co-linear, but instead are merely covariant?<jupyter_code>def plot_many(points):
    plt.scatter(*points.T)
    plt.xlim(0,25)
    plt.ylim(0,50)
    
points = np.random.multivariate_normal([15,15],[[12,6],[1,6]],50)
plot_many(points)
x, y = points.T
p2 = np.polyfit(x,y,1)
print(p2)
def f2(x):
    return p2[0] * x + p2[1]
plot_many(points)
nx = np.linspace(0,25,26)
plt.plot(nx,f2(nx),c='r');<jupyter_output><empty_output><jupyter_text>Regression by Gradient DescentLinear regression is implemented by a multitude of libraries:[A bunch of ways](https://medium.freecodecamp.org/data-science-with-python-8-ways-to-do-linear-regression-and-measure-their-speed-b5577d75f8b)You will never need to implement linear regression by hand, except right now! Linear regression is great because it has an analytic solution to the system of equations $Ax = b$, which can be computed via matrix multiplies:$x = (A^TA)^{-1}A^Tb$However, it is also trivial to compute the same solution via gradient descent. Remember those function gradients we studied at the beginning of class?<jupyter_code>points = np.random.multivariate_normal([10,10],[[10,10],[0,10]],40)
x,y = points.T
solution = np.polyfit(x,y,1)
plt.plot(x,y,marker='.',lineStyle="None");
plt.plot(x,x*solution[0] + solution[1],c='r');
plt.title('Cool line bro');<jupyter_output>/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: covariance is not positive-semidefinite.
  """Entry point for launching an IPython kernel.<jupyter_text>The solution is just a pair of scalars, $m$ and $b$. What if we pick them at random?<jupyter_code>m = np.random.random()
b = np.random.random()

# How is this solution?
plt.plot(x,y,marker='.',lineStyle="None");
plt.plot(x,x*solution[0] + solution[1],c='r');
plt.plot(x,m*x + b,c='g');<jupyter_output><empty_output><jupyter_text>Root-mean Squared Error (RMSE)RMSE is the most common "loss function", that is, a way of measuring the error in a model prediction. The error of a prediction is the difference between the values it predicted and the values that were expected.RMSE is calculated:$RMSE = \sqrt{\frac{\sum_{t=1}^T \left( \hat{y_t} - y_t \right) ^2}{T}}$$\hat{y_t}$ is the prediction of the model for a particular sample $t$. $y_t$ is the true label for sample $t$.Linear regression is the proess of minimizing the least-squared error between the line-of-regression and the sample points.If we calculate `polyfit` as above using `full=True`, it will return the error of the model computed by polyfit.<jupyter_code>coefficients, residual, solution_rank, solution_singular_values, rcond_epsilon = np.polyfit(x,y,1,full=True)
print(coefficients)
print(residual)
print(solution_rank)
print(solution_singular_values)
print(rcond_epsilon)<jupyter_output>[0.5176806  5.12747689]
[371.15112768]
2
[1.40663078 0.1462527 ]
8.881784197001252e-15<jupyter_text>Of particular interest here is the residual - this number is the MSE of the solution returned by `polyfit`. Can we reproduce it?What about the residual of the line we just guessed at?---<jupyter_code>confirm_residual = np.sum((x*coefficients[0] + coefficients[1] - y)**2)
print(confirm_residual)

random_residual = np.sum((x*m+b-y)**2)
print(random_residual)<jupyter_output>371.15112767623697
506.5687005771543<jupyter_text>Loss functionThe goal of gradient descent is to minimize the loss function - in this case, the loss function being RMSE of our estimated line. Without `polyfit`, we'd try to minimize the loss function all the way to 0, but `polyfit` has already calculated the best possible result! $MSE = 419$, $RMSE = \sqrt{419}$. Gradient Descent`polyfit` determined that the best possible solution for $f(x) = mx+b$ is `coefficients[0] * x + coefficients[1]`. For any other value of `m` and `b`, the error (or RMSE), will be larger. The error function, RMSE is $\sqrt(mx+b-y)^2$<jupyter_code>def mse(x,prediction_m,prediction_b,true_m,true_b):
    width = len(prediction_m)
    height = len(prediction_m[0])
    result = np.zeros(width*height).reshape(width,height)
    for row in range(width):
        for col in range(height):
            # calculate the RMSE of prediction_m/b, store in result
            m = prediction_m[row][col]
            b = prediction_b[row][col]
            result[row][col] = np.sum(((x*m+b)-(x*true_m+true_b))**2)
    return result

m_range = np.linspace(coefficients[0]*0.5, 2*coefficients[0], 15)
b_range = np.linspace(coefficients[1]*0.5, 2*coefficients[1], 15)
mm, bb = np.meshgrid(m_range,b_range)
Z = mse(x,mm,bb,coefficients[0],coefficients[1])

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D, get_test_data
from matplotlib import cm
import numpy as np

fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(mm, bb, Z,  cmap=cm.coolwarm,\
                       linewidth=0, antialiased=False)
ax.set_xlabel('m')
ax.set_ylabel('b')
ax.set_zlabel('residual')
ax.set_title('Gradient of Linear Regression');
print(coefficients)<jupyter_output>[0.5176806  5.12747689]<jupyter_text>Gradient Descent by Partial DerivativeThe error function $E_{rms} = \sum{\sqrt{(\hat y_i - y_i)^2}}$ is differentiable: that is, we can find the derivative of the error function, which informs us how to update the parameters to minimize it. A model with zero error is perfect.One of the foundations of machine learning is taking the partial derivative of the error function. When the error function is well defined like above, the gradient can be computed analytically. When it is not well defined, it can be computed using the method of finite differences which we studied in week 1. Lets drop the $\sqrt{}$ to make the derivative easier:$E_{rms} = \sum{(\hat y_i - y_i)^2} = \sum{(mx_i + b - y_i)^2}$We can't drop the $^2$, because that's what makes this a minimization problem.We have two variables $m$ and $b$ that we'd like to solve for:$\frac{\partial}{\partial m} = -2\sum x(mx_i + b - y)$$\frac{\partial}{\partial b} = -2\sum (mx_i + b - y)$See* [ML Cheatsheet](http://ml-cheatsheet.readthedocs.io/en/latest/linear_regression.html)and* [MCCormickML](http://mccormickml.com/2014/03/04/gradient-descent-derivation/)<jupyter_code>def linear_regression(X, y, m_current=0, b_current=0, epochs=10000, learning_rate=0.0001):
    N = float(len(y))
    for i in range(epochs):
        y_current = (m_current * X) + b_current
        cost = sum([data**2 for data in (y-y_current)]) / N
        m_gradient = -(2/N) * sum(X * (y - y_current))
        b_gradient = -(2/N) * sum(y - y_current)
        m_current = m_current - (learning_rate * m_gradient)
        b_current = b_current - (learning_rate * b_gradient)
    return m_current, b_current, cost

m_solution = linear_regression(x, y)
print(m_solution)
plt.plot(x,y,'o',marker='.')
plt.plot(x,x*solution[0]+solution[1],c='r')
plt.plot(x,x*m_solution[0]+m_solution[1],c='y');<jupyter_output>(0.9395646833921867, 0.4946374653798378, 10.187165644142938)<jupyter_text>Your AssignmentUse gradient descent on the Titanic dataset.* 1) Load and clean the data* 1.1) Pay special attention to outliers. How will your model perform with and without outliers?* 2) Identify which features have continuous numeric values* 3) Compute the linear regression of these values (there are only two predictable labels)* 3.1) Choose a set of features (columns) from Titanic and see which features minimize your RMSE* 4) Report on your RMSE (residual error) for each of the models that you produce.Linear regression works well with univariate data (that is, $dim(x) = 1 \times n$ and $dim(y) = 1 \times n$) or with multivariate data: $dim(x)) = m\times n$. Report your best regression selected from the data in the univariate case and the multivarate case. How do they differ?* 5) Would this model successfully produce good predictions on the chosen labels?  Thinking about your assignment1) What am I being asked to do?Perform feature selection where the selection criteria is the minimum RMSE that a linear regression algorithm will find. Do this for the univariate case - in some sense answering which feature has the most signal, assuming a linear relationship - and for the multivariate case.For this assignment, the labels I will try to predict will be the two continuous variables: `fare` and `age`. Note that this is different from the canonical label that models typically try to predict - `survived`.2) What coding steps need to be taken to satisfy the problem?The common data preparation steps - dropping duplicates, treating missing values, normalizing data if appropriate, followed by linear regression on the continuous variables - `fare` and `age` - and error calculation for all combinations of features.3) What must I do to claim that I have "completed" the assignment?* Relate subsets of Titanic features to the RMSE of a linear model for the two continuous labels* Identify the feature that produces the model with least RMSE out of all sets of single features* Identify the set of features that produce the model with the least RMSE out of all possible sets of features* Discuss the quality of the predictions for the best models found* Repeat the above analysis on two versions of the dataset, one without outliers and one with<jupyter_code># Complete assignment here.

import pandas as pd
import seaborn as sns
from sklearn.preprocessing import LabelEncoder, OneHotEncoder
from sklearn.linear_model import LinearRegression
from itertools import chain, combinations<jupyter_output><empty_output><jupyter_text>Load and Preprocess Data<jupyter_code>data = sns.load_dataset('titanic')
data = data.drop(['alive','adult_male','who','class','embark_town', 'deck'], axis=1)
data = data.drop_duplicates()

data.head()
def fill_mixed_median_mode(dataframe, medians=list()):
    """ Fill missing values with median for specified column, otherwise mode
    
    Args:
        dataframe (pandas.core.frame.DataFrame): rows of observations of features
        medians (list): columns to fill missing values with median instead of mode
        
    Returns:
        dataframe with no missing values
    """
    
    
    null = dataframe.isnull().any()
    null_cols = list(null[null].index)
    
    fill = pd.Series([data[c].median() if c in medians else data[c].mode()[0]
                     for c in null_cols], index=null_cols)
    
    dataframe[null_cols] = dataframe[null_cols].fillna(fill)
    return dataframe
data_f = fill_mixed_median_mode(data, ['Age'])
for label in ['embarked','sex', 'alone']:
    data_f[label] = LabelEncoder().fit_transform(data_f[label])

embarked_one_hot = OneHotEncoder().fit_transform(data_f[['embarked']]).toarray()
embarked = pd.DataFrame(embarked_one_hot, 
                        columns=['Southampton', 'Cherbourg', 'Queenstown'], 
                        dtype=np.int64)

data_f = data_f.reset_index(drop=True)
data_enc = data_f.join([embarked])
data_enc = data_enc.drop(['embarked'], axis=1)

data_enc.head()<jupyter_output><empty_output><jupyter_text>I will include `survived` in my set of features. My two outcomes are `age` and `fare`. Every other feature will also be considered.  Univariate Linear Regression<jupyter_code>def compare_models(data, features, outcomes):
    models = pd.DataFrame([])
    
    for feature in features:
        if type(feature) != list:
            feature = [feature]
        
        for outcome in outcomes:
            X = data[feature].values.reshape(-1, len(feature))
            Y = data[outcome].values.reshape(-1, 1)
            
            model = LinearRegression()
            model.fit(X, Y)
            yhat = model.predict(X)
            
            rmse = np.sqrt(np.sum(np.power(yhat-Y, 2))/yhat.shape[0])
            params = model.coef_[0]
            intercept, = model.intercept_
            
            summary = pd.Series([outcome, feature, rmse, params, intercept], name='')
            models = models.append(summary)
            
    models.columns = ['outcome', 'features', 'RMSE', 'parameters', 'intercept']
    
    return models
features = data_enc.drop(['age', 'fare'], axis=1).columns
outcomes = ['age', 'fare']

models_uni = compare_models(data_enc, features, outcomes)
models_uni = models_uni.sort_values(by='RMSE')

age_uni = models_uni[models_uni['outcome']=='age'].drop('outcome', axis=1)
fare_uni = models_uni[models_uni['outcome']=='fare'].drop('outcome', axis=1)

age_uni = age_uni.reset_index(drop=True)
fare_uni = fare_uni.reset_index(drop=True)
age_uni
fare_uni<jupyter_output><empty_output><jupyter_text>It seems that for both `age` and `fare`, `pclass` is the feature that produces the model with the least RMSE. However, one should keep in mind that this is on unnormalized data. A better metric for the quality of a linear model would be its $R^2$ score.  Multivariate Linear Regression<jupyter_code>def all_subsets(data_cols):
    return chain(*map(lambda x: combinations(data_cols, x), range(0, len(data_cols)+1)))
subsets = [list(subset) for subset in all_subsets(features) if len(subset) > 0]
subsets[:15]
models = compare_models(data_enc, subsets, outcomes)
age = models[models['outcome']=='age'].drop('outcome', axis=1)
fare = models[models['outcome']=='fare'].drop('outcome', axis=1)

age = age.reset_index(drop=True).sort_values(by='RMSE')
fare = fare.reset_index(drop=True).sort_values(by='RMSE')
age.head()
age['features'].iloc[0]
fare.head()
fare['features'].iloc[0]<jupyter_output><empty_output><jupyter_text>For both `fare` and `age`, a model that considered 8 features obtained the lowest RMSE. These are the features `survived`, `pclass`, `sex`,  `sibsp`, `parch` and `alone` along with two out of the three cities. Interestingly, the cities are not the same for predicting each outcome.In the future, I could consider including one less one-hot feature than was encoded, since its information is completely contained in the baseline case of all the other one-hot features being equal to zero. This would remedy the "dummy variable trap" contributing to [multicollinearity](https://en.wikipedia.org/wiki/MulticollinearityRemedies_for_multicollinearity).  Discussion<jupyter_code>fig, ax = plt.subplots(2)
fig.subplots_adjust(hspace=0.3)

for ix, outcome in enumerate(outcomes):
    line = models_uni[models_uni['outcome']==outcome].iloc[0]
    feature = line['features'][0]
    
    x = data_enc[feature]
    y = data_enc[outcome]
    ax[ix].scatter(x, y)
    x_line = np.linspace(x.min(), x.max())
    y_line = x_smooth*line['parameters'][0] + line['intercept']
    ax[ix].plot(x_line, y_line)
    
    ax[ix].set(xlabel=feature, ylabel=outcome);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1173485176703596, 'fraction_numerical': 0.023519217676548864, 'mean_word_length': 4.6774420238931835, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 57, 'https://': 2, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '404230', 'license': 'MIT', 'n_tokens_mistral': 5538, 'n_tokens_neox': 5171, 'n_words': 1663, 'path': 'Week 06 Modeling/Linear Regression.ipynb', 'repo_name': 'jraval/LambdaSchoolDataScience', 'size': '16157'}