<jupyter_start><jupyter_text>Python Practice Lecture 13 MATH 342W Queens College - Polynomial Regression Author: <NAME> Date: March 15, 2021 Agenda:* Nonlinear Regression: Polynomial Regression* Extrapolation vs Interpolation* Overfitting with Polynomials Nonlinear Linear Regression: Polynomial RegressionLet's generate a polynomial model of degree 2, that is, $f = h^* \in \mathcal{H}$ and let $\epsilon$ be random noise (the error due to ignorance) for $\mathbb{D}$.<jupyter_code># importing dependencies
import numpy as np
from sklearn.linear_model import LinearRegression
import pandas as pd
import seaborn as sns
sns.set(rc={"figure.figsize":(8, 5)}) #width=8, #height=5
sns.set_palette(sns.color_palette("colorblind")) # setting color palette
# Set seed
np.random.seed(1640)

# Generating h_*
n = 30
beta_0 = 1
beta_1 = 0
beta_2 = 1

x = np.random.uniform(low=-2.0, high=5.0, size=n)

# best possible model
h_star_x = beta_0 + beta_1 * x + beta_2 * x**2

# actual data differs due to information we don't have
epsilon = [np.random.normal(0, 1) for i in range(n)]
y = h_star_x + epsilon

df = pd.DataFrame({'x':x, 'y':y, 'h_star_x':h_star_x})

sns.scatterplot(data = df, x = x, y = y)<jupyter_output><empty_output><jupyter_text>Let's try to estimate with a line:<jupyter_code># importing dependencies
from typing import Sequence, Union
import numpy as np
import matplotlib.pyplot as plt

# function to plot via slope and intercept.
def abline(slope, intercept):
    """Plot a line from slope and intercept"""
    axes = plt.gca()
    x_vals = np.array(axes.get_xlim())
    y_vals = intercept + slope * x_vals
    plt.plot(x_vals, y_vals, '--')


# functions to plot polynomial expressions
Number = Union[int, float, complex]

def polyval(coefficients: Sequence[Number], x: Sequence[Number]) -> np.ndarray:
    # expand dimensions to allow broadcasting (constant time + inexpensive)
    # axis=-1 allows for arbitrarily shaped x
    x = np.expand_dims(x, axis=-1)
    powers = x ** np.arange(len(coefficients))
    return powers @ coefficients

def polyplot(coefficients: Sequence[Number], x: Sequence[Number]) -> None:
    y = polyval(coefficients, x)
    plt.plot(x, y)
model = LinearRegression(fit_intercept = False)

model.fit(pd.DataFrame(df.x), pd.DataFrame(df.y))

slope = model.coef_[0][0]
intercept = model.intercept_

# plot data points
sns.scatterplot(data = df, x = x, y = y)

# plot model regression line
abline(slope, intercept)

# plotting h_star_x
polyplot(np.array([beta_0, beta_1, beta_2]), np.linspace(-2, 5, 300))
# output R^2 of dotted line
model.score(pd.DataFrame(df.x), pd.DataFrame(df.y))<jupyter_output><empty_output><jupyter_text>The relationship is "underfit". $\mathcal{H}$ is not rich enough right now to express something close to $f(x)$. But it is better than the null model!Now let's do a polynomial regression of degree two. Let's do so manually:<jupyter_code># Set seed
np.random.seed(1640)

# defining X
X = [[1, x[i], x[i]**2] for i in range(len(x))]
X = pd.DataFrame(X)

X.head()
# solve for your b_vec
b = np.linalg.inv(X.transpose() @ X) @ X.transpose() @ y
print(b)
print([beta_0, beta_1, beta_2])<jupyter_output>0    0.342843
1   -0.219180
2    1.104987
dtype: float64
[1, 0, 1]<jupyter_text>And we'll plot it.<jupyter_code># plot data points
sns.scatterplot(data = df, x = x, y = y)

# plot model regression line
#abline(slope, intercept)

# plotting h_star_x - Blue
polyplot(np.array([beta_0, beta_1, beta_2]), np.linspace(-2, 5, 300))

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(-2, 5, 300))<jupyter_output><empty_output><jupyter_text>Now let's try polynomial of degree 3:<jupyter_code># Set seed
np.random.seed(1640)

# defining X
X = [[1, x[i], x[i]**2, x[i]**3] for i in range(len(x))]
X = pd.DataFrame(X)

X.head()
# solve for your b_vec
b = np.linalg.inv(X.transpose() @ X) @ X.transpose() @ y
print(b)
# plot data points
sns.scatterplot(data = df, x = x, y = y)

# plot model regression line
#abline(slope, intercept)

# plotting h_star_x - Blue
polyplot(np.array([beta_0, beta_1, beta_2]), np.linspace(-2, 5, 300))

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(-2, 5, 300))<jupyter_output><empty_output><jupyter_text>Still the same. Why? The $x^3$ term is like adding one "nonsense" predictor. One nonsense predictor marginally affects $R^2$ but it doesn't matter too much.Now let's try polynomial of degree 8:<jupyter_code># Set seed
np.random.seed(1640)

# defining X
X = [[1, x[i], x[i]**2, x[i]**3,  x[i]**4, x[i]**5, x[i]**6, x[i]**7, x[i]**8] for i in range(len(x))]
X = pd.DataFrame(X)

X.head()
# solve for your b_vec
b = np.linalg.inv(X.transpose() @ X) @ X.transpose() @ y
print(b)
# plot data points
sns.scatterplot(data = df, x = x, y = y)

# plot model regression line
#abline(slope, intercept)

# plotting h_star_x - Blue
polyplot(np.array([beta_0, beta_1, beta_2]), np.linspace(-2, 5, 300))

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(-2, 5, 300))<jupyter_output><empty_output><jupyter_text>We are seeing now a little bit of "overfitting" in the edge(s). We now have $p=9$ and $n=100$. We can do a lot worse!Now let's do polynomial of degree 13.<jupyter_code># Set seed
np.random.seed(1640)

# defining X
X = [[1, x[i], x[i]**2, x[i]**3,  x[i]**4, x[i]**5, x[i]**6, x[i]**7, x[i]**8,
      x[i]**9, x[i]**10, x[i]**11, x[i]**12, x[i]**13] for i in range(len(x))]

X = pd.DataFrame(X)

X.head()
# solve for your b_vec
b = np.linalg.inv(X.transpose() @ X) @ X.transpose() @ y
print(b)
# plot data points
sns.scatterplot(data = df, x = x, y = y)

# plot model regression line
#abline(slope, intercept)

# plotting h_star_x - Blue
polyplot(np.array([beta_0, beta_1, beta_2]), np.linspace(-2, 6, 300))

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(-2, 5, 300))<jupyter_output><empty_output><jupyter_text>What's happening for small values of $x$ (and a bit for large values)? This is called [Runge's Phenomenon](https://en.wikipedia.org/wiki/Runge%27s_phenomenon) meaning that the boundary activity of high-order polynomials has very large derivatives. Extrapolation vs InterpolationLet's take a look at the Galton Data used by Galton in 1885 to study the relationship between a parent's height and their childrens. We'll import it straight from a URL using Pandas' `.read_csv()` method.<jupyter_code># URL for Galton Height Data
url = "https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/HistData/Galton.csv"
galton_df = pd.read_csv(url)

galton_df = galton_df.iloc[:, 1:3]

galton_df.head()
model_g = LinearRegression(fit_intercept = True)

x = pd.DataFrame(galton_df["parent"])
y = pd.DataFrame(galton_df["child"])

model_g.fit(x, y)

b_0 = model_g.intercept_
b_1 = model_g.coef_[0]

print(b_0, b_1)
sns.scatterplot(data = galton_df, x = "parent", y = "child")

# plot model regression line
abline(b_1, b_0)<jupyter_output><empty_output><jupyter_text>The `Seaborn` jitter functionality is currently unavailable for scatter plots, so we'll create our own function for it.<jupyter_code># add jitter to object
def jitter(values,j = 0):
    return values + np.random.normal(j, 0.1, values.shape)

sns.scatterplot(x = jitter(galton_df["parent"]), y = jitter(galton_df["child"]),
                alpha = 0.5,  color='r')

# plot model regression line
abline(b_1, b_0)<jupyter_output><empty_output><jupyter_text>Let's say I want to predict child's height for parents' average height of 70in. All I do is:<jupyter_code>model_g.predict(pd.DataFrame({'parent':[70]}))[0][0]<jupyter_output><empty_output><jupyter_text>What if I want to predict for a parents' height of 5in. Is there any 12in tall human being? No... it is absurd. But nothing stops you from doing:<jupyter_code>model_g.predict(pd.DataFrame({'parent':[5]}))[0][0]<jupyter_output><empty_output><jupyter_text>That's [actually possible](https://www.guinnessworldrecords.com/news/2012/2/shortest-man-world-record-its-official!-chandra-bahadur-dangi-is-smallest-adult-of-all-time/).Look at our linear model from Euclid's perspective:<jupyter_code># plot data points
sns.scatterplot(x = jitter(galton_df["parent"]), y = jitter(galton_df["child"]),
                alpha = 0.4,  color='r')

b = [model_g.intercept_[0], model_g.coef_[0][0]]

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(-20, 100, 300))<jupyter_output><empty_output><jupyter_text>What is a linear model with $p + 1 = 2$. It's just a line. When geometry was first formalised by Euclid in the Elements, he defined a line to have "breadthless length" with a straight line being a line "which lies evenly with the points on itself". By "breadthless" he meant infinite in either direction. There is no mathematical problem with predicting childrens' heights using negative parents' heights e.g.<jupyter_code>model_g.predict(pd.DataFrame({'parent':[-5]}))[0][0]<jupyter_output><empty_output><jupyter_text>But this is absurd. So now we need to talk about a fundamental concept in data science we've been kind of avoiding and one that most people ignore. There are two types of prediction: interpolation and extrapolation. Interpolation is essentially the type of "prediction" we've been talking about this whole class. Extrapolation is totally different. It's what happens when you predict outside of the range of the covariate data you've seen in $\mathbb{D}$. Extrapolation is very dangerous - your models only work based on $\mathbb{D}$. Extrapolation is prediction outside of the range you've seen before which means you better have a good theoretical reason as to why your $\mathcal{H}$ function class will extend outside that range. Because each $\mathcal{H}$ function class will extrapolate very very differently. Overfitting with PolynomialsCan we achieve $R^2 = 100\%$ using polynomial regression? Yes. Here's an example in one dimension. These are called "interpolation polynomials". In one dimension, as long as the $x$ values are distinct, $n$ data points can be fit by a $n - 1$ degree polynomial. Here's an example with the same dataset we used in the beginning of the demo but with $n = 5$ instead.<jupyter_code># Set seed
np.random.seed(1921)

# Generating h_*
n = 5
beta_0 = 1
beta_1 = 0
beta_2 = 1

x = np.random.uniform(low=0, high=1, size=n)

# best possible model
h_star_x = beta_0 + beta_1 * x + beta_2 * x**2

# actual data differs due to information we don't have
epsilon = [np.random.normal(0, 1) for i in range(n)]
y = h_star_x + epsilon

df = pd.DataFrame({'x':x, 'y':y, 'h_star_x':h_star_x})

# plot it
sns.scatterplot(data = df, x = x, y = y)<jupyter_output><empty_output><jupyter_text>Now fit the polynomial model.<jupyter_code># Set seed
np.random.seed(1921)

# defining X
X = [[1, x[i], x[i]**2, x[i]**3, x[i]**4] for i in range(len(x))]
X = pd.DataFrame(X)

# solve for b_vec
b = np.linalg.inv(X.transpose() @ X) @ X.transpose() @ y

# plot data points
sns.scatterplot(data = df, x = x, y = y)

# plotting model polynomial - Orange
polyplot(np.array(b), np.linspace(0.24, 0.73, 100))<jupyter_output><empty_output><jupyter_text>Perfect fit! Observe the arc on the polynomial model above. All the values of $y$ are Normal mean 0, variance 1. If one provides an $x_* = 0.33$ the model would approximately return -19, which is far from a reasonable realization of our $\mathcal{N}(0, 1)$. Runge's Phenomenon has powerful effects.<jupyter_code># R^2
model = LinearRegression(fit_intercept = False)
model.fit(X, y)
model.score(X, y)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13399374348279458, 'fraction_numerical': 0.026068821689259645, 'mean_word_length': 5.031970649895178, 'pattern_counts': {'":': 1, '<': 65, '<?xml version=': 0, '>': 67, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, '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': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '383963', 'license': 'MIT', 'n_tokens_mistral': 4370, 'n_tokens_neox': 4107, 'n_words': 1280, 'path': 'Practice Notes/13 - Python Practice Lecture 13 MATH 342W Queens College - Polynomial Regression.ipynb', 'repo_name': 'Hi-Tree/342-2022-Python-Practice-lectures', 'size': '11515'}