<jupyter_start><jupyter_text>California housing dataset with linear and polynomial regression In this notebook, we'll use [linear regression](https://scikit-learn.org/stable/modules/linear_model.htmlordinary-least-squares), [regularized linear regression](https://scikit-learn.org/stable/modules/linear_model.htmlridge-regression), and [polynomial regression](https://scikit-learn.org/stable/modules/linear_model.htmlpolynomial-regression-extending-linear-models-with-basis-functions) to estimate median house values on Californian housing districts.First, the needed imports.<jupyter_code>%matplotlib inline

import numpy as np
from sklearn import datasets, __version__
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()<jupyter_output><empty_output><jupyter_text>DataThen we load the California housing data. First time we need to download the data, which can take a while.<jupyter_code>chd = datasets.fetch_california_housing()<jupyter_output><empty_output><jupyter_text>Let's first convert the data into a Pandas DataFrame to inspect some basic statistics:<jupyter_code>df = pd.DataFrame(data=chd.data, columns=chd.feature_names)
df['Target'] = pd.Series(chd.target, index=df.index)
df.describe()<jupyter_output><empty_output><jupyter_text>We see that the data consists of 20640 housing districts, each characterized with 8 attributes: *MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude*. We also defined  a target value (median house value) for each housing district. Let's then plot all attributes against the target value:<jupyter_code>plt.figure(figsize=(15,10))
for i in range(8):
    plt.subplot(4,2,i+1)
    plt.scatter(chd.data[:,i], chd.target, s=2, label=chd.feature_names[i])
    plt.legend(loc='best')<jupyter_output><empty_output><jupyter_text>We'll now split the data into a training and a test set. Let's use 5000 samples as test data.Let's also select a single attribute to start the analysis with, for example *MedInc*. This way we can plot the regression functions against the target value. Later we will use all attributes in the regression.<jupyter_code>test_size = 5000
single_attribute = 'MedInc'

X_train_all, X_test_all, y_train, y_test = train_test_split(
    chd.data, chd.target, test_size=test_size, shuffle=True)

attribute_index = chd.feature_names.index(single_attribute)
X_train_single = X_train_all[:, attribute_index].reshape(-1, 1)
X_test_single = X_test_all[:, attribute_index].reshape(-1, 1)
     
print()
print('California housing data: train:',len(X_train_all),'test:',len(X_test_all))
print()
print('X_train_all:', X_train_all.shape)
print('X_train_single:', X_train_single.shape)
print('y_train:', y_train.shape)
print()
print('X_test_all', X_test_all.shape)
print('X_test_single', X_test_single.shape)
print('y_test', y_test.shape)<jupyter_output><empty_output><jupyter_text>The training data matrix `X_train_all` is a matrix of size (`n_train`, 8), and `X_train_single` contains only the first attribute (*MedInc* by default) of each housing district. The vector `y_train` contains the target value (median house value) for each housing district in the training set.Let's start our analysis with the single attribute. Later, you can set `only_single_attribute = False` to use all eight attributes in the regression. As the final step, let's scale our input data to zero mean and unit variance.<jupyter_code>only_single_attribute = True

if only_single_attribute:
    X_train = X_train_single
    X_test = X_test_single
else:
    X_train = X_train_all
    X_test = X_test_all

scaler = StandardScaler().fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)
print('X_train: shape:', X_train.shape, 'mean:', X_train.mean(axis=0), 'std:', X_train.std(axis=0))
print('X_test: shape:', X_test.shape, 'mean:', X_test.mean(axis=0), 'std:', X_test.std(axis=0))<jupyter_output><empty_output><jupyter_text>Linear regressionWe begin with linear regression:$$J(w) = \|y - Xw\|^2_2$$ LearningThe parameters of linear regression can be solved in closed form as:$$\hat{w} = (X^TX)^{-1}X^Ty$$<jupyter_code>%%time

lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)
print('coefficients:', lin_reg.coef_)
print('intercept:', lin_reg.intercept_)<jupyter_output><empty_output><jupyter_text>We can visualize the results if we are using only a single attribute:<jupyter_code>if X_train.shape[1] == 1:
    plt.figure(figsize=(10, 10))
    plt.scatter(X_train, y_train, s=5)
    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)
    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear')
    plt.legend(loc='best');<jupyter_output><empty_output><jupyter_text>InferenceWe use *mean squared error* as the performance measure for our regression algorihm:<jupyter_code>%%time

predictions = lin_reg.predict(X_test)
print("Mean squared error: %.3f"
      % mean_squared_error(y_test, predictions))<jupyter_output><empty_output><jupyter_text>Regularized linear regression: RidgeRidge regression adds $L_2$ regularization: $$J(w) = \|y - Xw\|^2_2 + \alpha \|w\|^2_2$$where $\alpha \ge 0$ is the penalty parameter for the weights. You can experiment with different values of $\alpha$.You can also try out [`Lasso()`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.htmlsklearn-linear-model-lasso) or [`ElasticNet()`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.htmlsklearn-linear-model-elasticnet). Note that Elastic net has also a second parameter `l1_ratio`.  Learning<jupyter_code>%%time

alpha = 1e3

rdg_reg = Ridge(alpha=alpha)
rdg_reg.fit(X_train, y_train)
print('coefficients:', rdg_reg.coef_)
print('intercept:', rdg_reg.intercept_)
if X_train.shape[1] == 1:
    plt.figure(figsize=(10, 10))
    plt.scatter(X_train, y_train, s=5)
    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)
    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear');
    plt.scatter(reg_x, rdg_reg.predict(reg_x), s=8, label=r'ridge, $\alpha=${}'.format(alpha))
    plt.legend(loc='best');<jupyter_output><empty_output><jupyter_text>Inference<jupyter_code>%%time

predictions = rdg_reg.predict(X_test)
print("Mean squared error: %.3f"
      % mean_squared_error(y_test, predictions))<jupyter_output><empty_output><jupyter_text>Polynomial regressionPolynomial regression can be performed by constructing polynomial features, e.g.:$$z=[1,\,x_1,\,x_2,\,x_1x_2,\,x_1^2,\,x_2^2]$$and using a linear model with the new features:$$J(w) = \|z - X'w\|^2_2$$ LearningTo implement polynomial regression, we use scikit-learn's [Pipeline](https://scikit-learn.org/stable/modules/compose.htmlpipeline), a tool for building composite estimators. Note that the polynomial features contain all possible combinations, so the number of features grows quickly especially when using many attributes. Also, you can try using regularized linear regression with polynomial features.<jupyter_code>%%time
degree = 5

poly_model = Pipeline([('poly', PolynomialFeatures(degree=degree)),
                      ('linear', LinearRegression(fit_intercept=False))])
poly_model.fit(X_train, y_train)
print('coefficients:', poly_model.steps[1][1].coef_)
print('intercept:', poly_model.steps[1][1].intercept_)
if X_train.shape[1] == 1:
    plt.figure(figsize=(10, 10))
    plt.scatter(X_train, y_train, s=5)
    reg_x = np.arange(np.min(X_train), np.max(X_train), 0.01).reshape(-1, 1)
    plt.scatter(reg_x, lin_reg.predict(reg_x), s=8, label='linear');
    plt.scatter(reg_x, poly_model.predict(reg_x), s=8, label='polynomial')
    plt.legend(loc='best');<jupyter_output><empty_output><jupyter_text>Inference<jupyter_code>%%time

predictions = poly_model.predict(X_test)
print("Mean squared error: %.3f"
      % mean_squared_error(y_test, predictions))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12579617834394904, 'fraction_numerical': 0.012738853503184714, 'mean_word_length': 6.574211502782932, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '57675', 'license': 'MIT', 'n_tokens_mistral': 2878, 'n_tokens_neox': 2750, 'n_words': 652, 'path': 'notebooks/sklearn-chd-lr.ipynb', 'repo_name': 'CSCfi/machine-learning-scripts', 'size': '8164'}