<jupyter_start><jupyter_text>Intro to SML This is based on the notebook file [01 in <NAME>'s github page](https://github.com/ageron/handson-ml)<jupyter_code># Import the necessary packages
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')

# Where to save the figures
PROJECT_ROOT_DIR = "../.."
datapath = PROJECT_ROOT_DIR + "/data/lifesat/"<jupyter_output><empty_output><jupyter_text>Load and prepare Life satisfaction data<jupyter_code># Download CSV from http://stats.oecd.org/index.aspx?DataSetCode=BLI
oecd_bli = pd.read_csv(datapath+"oecd_bli_2015.csv", thousands=',')
oecd_bli = oecd_bli[oecd_bli["INEQUALITY"]=="TOT"]
oecd_bli.head()
oecd_bli = oecd_bli.pivot(index="Country", columns="Indicator", values="Value")
oecd_bli.head()
oecd_bli.columns
oecd_bli["Life satisfaction"].head()<jupyter_output><empty_output><jupyter_text>Load and prepare GDP per capita data<jupyter_code># Download data from http://goo.gl/j1MSKe (=> imf.org)
gdp_per_capita = pd.read_csv(datapath+"gdp_per_capita.csv", thousands=',', delimiter='\t',
                             encoding='latin1', na_values="n/a")
gdp_per_capita.rename(columns={"2015": "GDP per capita"}, inplace=True)
gdp_per_capita.set_index("Country", inplace=True)
gdp_per_capita.head()
full_country_stats = pd.merge(left=oecd_bli, right=gdp_per_capita, left_index=True, right_index=True)
full_country_stats.sort_values(by="GDP per capita", inplace="True")
full_country_stats.head()
full_country_stats[["GDP per capita", 'Life satisfaction']].loc["United States"]
remove_indices = [0, 1, 6, 8, 33, 34, 35]
keep_indices = list(set(range(36)) - set(remove_indices))

sample_data = full_country_stats[["GDP per capita", 'Life satisfaction']].iloc[keep_indices]
missing_data = full_country_stats[["GDP per capita", 'Life satisfaction']].iloc[remove_indices]
_ = sample_data.plot("GDP per capita",'Life satisfaction',kind='scatter')<jupyter_output><empty_output><jupyter_text>Linear regression to predict Life satisfaction using GDP per capita Let's make a design matrix and response vector.<jupyter_code>x = np.array(sample_data[['GDP per capita']])
n,_ = x.shape
X = np.hstack((np.ones((n,1)),x))
y = np.array(sample_data[['Life satisfaction']])
print(X[0:5,:])<jupyter_output>[[  1.00000000e+00   9.05491400e+03]
 [  1.00000000e+00   9.43737200e+03]
 [  1.00000000e+00   1.22398940e+04]
 [  1.00000000e+00   1.24953340e+04]
 [  1.00000000e+00   1.59917360e+04]]<jupyter_text>Ordinary least squares is defined by the normal equation$$X^\top X \hat \beta = X^\top y$$$$\hat y = X \hat \beta$$Let's fit the linear model and predict the points.<jupyter_code>from sklearn import linear_model
lin1 = linear_model.LinearRegression(fit_intercept=False)
lin1.fit(X,y)
print(lin1.coef_)
yhat = lin1.predict(X)
# Compute the training error
RSS = np.sum((y - yhat)**2.)
training_error = RSS / n
print(training_error)
_ = plt.scatter(X[:,1],y)
_ = plt.plot(X[:,1],yhat,'r')<jupyter_output><empty_output><jupyter_text>The regression line seems to demonstrate a good fit to a linear model.  We can use it to predict a new point, to assess the impact of a tax policy for example.  Quadratic Regression<jupyter_code>X_quad = np.hstack((np.ones((n,1)),x,x**2))<jupyter_output><empty_output><jupyter_text>In the above, we add a quadratic term to the design matrix, to see if it improves the fit.<jupyter_code>from sklearn import linear_model
lin1 = linear_model.LinearRegression(fit_intercept=False)
lin1.fit(X_quad,y)
yhat = lin1.predict(X_quad)
# Compute the training error
RSS = np.sum((y - yhat)**2.)
training_error = RSS / n
print(training_error)
_ = plt.scatter(X[:,1],y)
_ = plt.plot(X[:,1],yhat,'r')<jupyter_output><empty_output><jupyter_text>The fits are essentially the same between linear and quadratic fit, and the training risk is nearly equal.  Leave-one-out prediction simulation In this section, I show that the training risk is an optimistic measure of the true risk.<jupyter_code>from sklearn.model_selection import LeaveOneOut
loo = LeaveOneOut()
loo.get_n_splits(X)
loo_losses = []
for train_index, test_index in loo.split(X):
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]
    lin1.fit(X_train,y_train)
    y_hat = lin1.predict(X_test)
    loss = np.sum((y_hat - y_test)**2)
    loo_losses.append(loss)
    print("TEST: " + sample_data.index.values[test_index][0] + " " + str(loss))
plt.hist(loo_losses)
plt.vlines(training_error,0,12,'b',label="Training Risk")
plt.vlines(np.mean(loo_losses),0,12,'k',label="LOO Ave Risk")
_ = plt.legend()<jupyter_output><empty_output><jupyter_text>This histogram of the leave-one-out errors demonstrates that the training risk can be optimistic, because the Leave-one-out risk (the average of the LOO errors) is larger than the training risk.<jupyter_code>def loo_risk(X,y):
    """
    Construct the leave-one-out square error risk for the linear model
    
    Input: design matrix, X, response vector, y
    Output: scalar LOO risk
    """
    loo = LeaveOneOut()
    loo.get_n_splits(X)
    loo_losses = []
    for train_index, test_index in loo.split(X):
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]
        lin1.fit(X_train,y_train)
        y_hat = lin1.predict(X_test)
        loss = np.sum((y_hat - y_test)**2)
        loo_losses.append(loss)
    return np.mean(loo_losses)
help(loo_risk)
print("Linear LOO risk: {}".format(loo_risk(X,y)))
print("Quadratic LOO risk: {}".format(loo_risk(X_quad,y)))<jupyter_output>Linear LOO risk: 0.209592810168
Quadratic LOO risk: 0.230390450851
{'fraction_non_alphanumeric': 0.13503971756398941, 'fraction_numerical': 0.033715798764342456, 'mean_word_length': 5.697399527186761, 'pattern_counts': {'":': 1, '<': 36, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '66192', 'license': 'MIT', 'n_tokens_mistral': 2240, 'n_tokens_neox': 2088, 'n_words': 459, 'path': '2018_material/lectures/lecture1/lecture1.ipynb', 'repo_name': 'jynwang/DavisSML', 'size': '5673'}