PolynomialRegression/Polynomial Regression.ipynb

Algorytm najszybszego spadku dla regresji wielomianowej.

Skład grupy:

• Nowak Anna,
• Łaźna Patrycja,
• Bregier Damian

0. Podstawowe informacje o zbiorze danych

Ze względu na specyfikę wytycznych projektu zakładającego wykorzystanie jedynie dwóch cech: x i y oraz modelowanie zależności y od x za pomocą funkcji wielomianowej dobrany został przez nas odpowiedni dataset. Obejmuje on jedynie trzy kolumny, czyli: płeć, wzrost w calach oraz wagę w funtach, z czego ze względu na specyfikę projektu wykorzystywane są jedynie wzrost oraz waga. Każdy z parametrów zawiera po 10 tysięcy unikalnych wartości.

Dokładne informacje na temat tego zbioru, jego zawartości i samych danych można znaleźć pod adresem: https://www.kaggle.com/mustafaali96/weight-height.

import pandas as pd
import random
import matplotlib.pyplot as plt
import numpy as np
import ipywidgets as widgets
np.set_printoptions(suppress=True)

1. Wczytywanie i preprocessing danych

# Stopień wielomianu
degree = 4
plot_degree = 1
X_plot_lib = np.linspace(0, 100, 1000)
initial_theta = np.matrix([0] * (degree + 1)).reshape(degree + 1, 1)

# Wybór dwóch kolumn - dotyczących wzrostu i wagi
data = pd.read_csv('weight-height.csv')[["Height", "Weight"]]
# Czyszczenie tabeli i wartości pustych
data = data.dropna()
data_matrix = np.matrix(data)

m, n_plus_1 = data_matrix.shape
n = n_plus_1 - 1
X = (np.ones((m, 1)))

for i in range(1, degree + 1):
    Xn = np.power(data_matrix[:, 0:n], i)
    Xn /= np.amax(Xn, axis=0)
    X = np.concatenate((X, Xn), axis=1)

X = np.matrix(X).reshape(m, degree * n + 1)
Y = np.matrix(data_matrix[:, -1])

2. Metody regresji wielomianowej

Regresja wielomianowa służy do obliczania zależności pomiędzy zmienną zależną a jedną lub więcej zmiennymi niezależnymi, które mogą występować w wyższych potęgach. Wykorzystywana w przypadku zbiorów cechujacych się nieliniowymi zależnościami.

Równanie regresji wielomianowej:

def print_theta(theta):
    print("f(x) = ", end="")
    for i,x in enumerate(theta.tolist()[:-1]):
        x = x[0]
        print(f"{x}x^{i}", end=" + ")
    print(f"{theta.tolist()[-1][0]}^{len(theta) - 1}")

# Implementacja wzoru na regresję wielomianową
def polynomial_regression(theta, x):
    x = x/data["Height"].max()
    return sum(theta * np.power(x, i) for i, theta in enumerate(theta.tolist()))

# Implementacja wzoru na RMSE, czyli pierwiastek z błędu średniokwadratowego
def mean_squared_error(theta, X, Y):
    J = 1.0 / (2.0 * m) * ((X * theta - Y).T * (X * theta - Y))
    return J.item()

3. Wybrane rodzaje gradientów

Metoda gradientu prostego to algorytm numeryczny mający na celu znalezienie minimum lokalnego dla zadanej funkcji celu. Jak sama nazwa wskazuje jest to jedna z prostszych metod, która opiera się na iteracyjnym poszukiwaniu minimum. Proces ten (dla pierwszych czterech kroków), został pokazany na poniższym rysunku.

Metoda gradientu prostego

def gradient(theta, X, Y):
    return 1.0 / len(Y) * (X.T * (X * theta - Y)) 

Metoda najszybszego spadku jest iteracyjnym algorytmem wyszukiwania minimym zadanej funkcji celu i stanowi modyfikację gradientu prostego. Na samym początku algorytmu wybierany jest punkt startowy, w którym jest obliczany antygradient funkcji celu stanowiący kierunek poszukiwań w procedurze. Kolejny krok jest obliczany w oparciu o poprzedni i jeżeli nie spełnia on warunku stopu to obliczenia są powtarzane do skutku.

Podstawową różnicę w stosunku do gradientu prostego stanowi sposób wyszukiwania wartości, która jest dokonywana przez minimalizację kierunkową funkcji.

# Steepest descent
def steepest_descent(X, Y, theta, cost_function = mean_squared_error, 
                     single_step=0.01, epochs=30):
    cost = cost_function(theta, X, Y)
    logs = [[cost, theta]]
    for i in range(epochs):
        direction = gradient(theta, X, Y)
        j = 0
        while(True):
            j+=1
            next_theta = theta - single_step * direction
            next_cost = cost_function(next_theta, X, Y)
            if(next_cost >= cost or type(next_cost) != float):
                break
            cost = next_cost
            theta = next_theta
        logs.append([next_cost, theta])
    return theta, logs

Ponadto w toku prac projektowych przetestowane zostaly także inne algorytmy takie jak batch gradient descend (BGD), BGD z momentum, mini-batch gradient descend oraz stochastic gradient descend.

BGD to metoda w której błąd obliczany jest dla każdego parametru ze zbioru treningowego, ale sam model jest aktualizowany dopiero po ocenie wszystkich danych.

MBGD opiera się na podziale danch treningowych na niewielkie pakiety, które są potem wykorzystywane do obliczania wartości błędu i aktualizacji modelu.

# Batch gradient descend (BGD)
def BGD(X, Y, theta, cost_function = mean_squared_error, learning_rate=0.1, 
        eps=10**-5, max_steps = 10000000000):
    cost = cost_function(theta, X, Y)
    logs = [[cost, theta]]
    
    for i in range(max_steps):
        theta = theta - learning_rate * gradient(theta, X, Y)
        next_cost = cost_function(theta, X, Y)
        logs.append([next_cost, theta])
        if abs(cost - next_cost) <= eps:
            break
        cost = next_cost
    return theta, logs

# BGD z momentum
def momentum(X, Y, theta, cost_function = mean_squared_error, 
             learning_rate=0.1, momentum = 0.3, epochs=30, degree=degree):
    cost = cost_function(theta, X, Y)
    logs = [[cost, theta]]
    delta_history = [np.matrix([0]*(degree+1)).T]
    for i in range(epochs):
        delta_history.append(momentum * delta_history[-1] + 
                             learning_rate * gradient(theta, X, Y))
        theta = theta - delta_history[-1]
        next_cost = cost_function(theta, X, Y)
        logs.append([next_cost, theta])
        cost = next_cost
    return theta, logs

# Mini-batch gradient descend (MBGD)
def MBGD(X, Y, theta, cost_function = mean_squared_error, 
         learning_rate=0.1, epochs=5, batch_size=16):
    cost = cost_function(theta, X, Y)
    logs = [[cost, theta]]
    start, end = 0, batch_size
    
    steps = m / batch_size
    for i in range(epochs):
        zipped_XY = list(zip(X, Y))
        random.shuffle(zipped_XY)
        X_shuffled, Y_shuffled = zip(*zipped_XY)
        X_shuffled = np.concatenate(X_shuffled, axis=0) 
        Y_shuffled = np.concatenate(Y_shuffled, axis=0) 
        for j in range(int(steps)):
            batch =  X_shuffled[start:end,:], Y_shuffled[start:end,:]
            theta = theta - learning_rate * gradient(theta, batch[0], batch[1])
            cost = cost_function(theta, X, Y)
            logs.append([cost, theta])

            if start + batch_size < batch_size:
                start += batch_size
            else:
                start = 0
            end = min(start + batch_size, m)
    return theta, logs

# Stochastic gradient descend (SGD)
def SGD(X, Y, theta, cost_function = mean_squared_error, learning_rate=0.1,
        epochs=5):
    return MBGD(X, Y, theta, cost_function, learning_rate, epochs, 1)

4. Zestawienie wyników

final_theta_SGD, logs_3 = SGD(X, Y, initial_theta, epochs = 60)
print("SGD:")
print_theta(final_theta_SGD)
print(f"{len(logs_3)} updates\n")

final_theta_BGD, logs_1 = BGD(X, Y, initial_theta)
print("BGD:")
print_theta(final_theta_BGD)
print(f"{len(logs_1)} updates\n")

final_theta_MBGD, logs_2 = MBGD(X, Y, initial_theta, epochs = 60)
print("MBGD:")
print_theta(final_theta_MBGD)
print(f"{len(logs_2)} updates\n")


final_theta_momentum, logs_4 = momentum(X, Y, initial_theta, epochs = 60)
print("momentum:")
print_theta(final_theta_momentum)
print(f"{len(logs_4)} updates\n")

final_steepest_descent, logs_5 = steepest_descent(X, Y, initial_theta, epochs = 60)
print("steepest descent:")
print_theta(final_steepest_descent)
print(f"{len(logs_5)} updates\n")

SGD:
f(x) = 28.344374625703598x^0 + 39.622015878252434x^1 + 46.23584544945583x^2 + 49.503407925541254x^3 + 50.40885417329693^4
600001 updates

BGD:
f(x) = -65.64506277392525x^0 + 68.77648241256007x^1 + 114.75349304333506x^2 + 102.20675310280686x^3 + 52.800146520466946^4
82705 updates

MBGD:
f(x) = -46.04759399843719x^0 + 42.66139616280628x^1 + 87.3646802267618x^2 + 102.16290156608075x^3 + 97.32102412195839^4
37501 updates

momentum:
f(x) = 46.73334506478951x^0 + 46.041369573952764x^1 + 44.48349874963423x^2 + 42.37494564787237x^3 + 39.9436244821308^4
61 updates

steepest descent:
f(x) = 13.638648831593002x^0 + 39.09249361982211x^1 + 55.995416451700564x^2 + 66.58927282539214x^3 + 72.56809724003132^4
61 updates

Wniosek: Metoda najszybszego spadku oraz momentum cechują się najszybszą zbieżnością do minimum lokalnego danej funkcji. Przeciętnie, dla zadanych parametrów, minimum jest osiągane w ciągu 61 aktualizacji. Podczas gdy SGD wymaga ich ponad 600 tysięcy.

5. Reprezentacja graficzna regresji wielomianowej i funkcji kosztu

def plot_polynomial_regression(theta):
    fig = plt.figure(figsize=(10,5))
    Y_plot = polynomial_regression(theta, X_plot_lib).tolist()
    chart = fig.add_subplot()
    chart.plot(data["Height"], Y ,"go")
    chart.plot(X_plot_lib, Y_plot, color="red", lw=2, label=f"degree {len(theta)}")
    plt.ylim([0,250])
    plt.show()
    

def plot_cost_function(logs):
    no_iter = no_iter = [x for x in range(1, len(logs) + 1)]
    all_cost = [row[0] for row in logs]
    
    fig = plt.figure(figsize=(10,5))  
    plt.plot(no_iter, all_cost)
    plt.xlabel('No. of iterations')
    plt.ylabel('Cost')
    plt.show()
    

print("BGD:")
plot_polynomial_regression(final_theta_BGD)
plot_cost_function(logs_1)
print("MBGD:")
plot_polynomial_regression(final_theta_MBGD)
plot_cost_function(logs_2)
print("SGD:")
plot_polynomial_regression(final_theta_SGD)
plot_cost_function(logs_3)
print("Momentum:")
plot_polynomial_regression(final_theta_momentum)
plot_cost_function(logs_4)
print("Metoda najszybszego spadku:")
plot_polynomial_regression(final_steepest_descent)
plot_cost_function(logs_5)

BGD:

MBGD:

SGD:

Momentum:

Metoda najszybszego spadku:

Wizualizacja metody najszybszego spadku dla wielomianu drugiego stopnia

initial_theta_plot = np.matrix([0] * (plot_degree + 1)).reshape(plot_degree + 1, 1)
m, n_plus_1 = data_matrix.shape
n = n_plus_1 - 1
X_plot = (np.ones((m, 1)))

for i in range(1, plot_degree + 1):
    Xn = np.power(data_matrix[:, 0:n], i)
    Xn /= np.amax(Xn, axis=0)
    X_plot = np.concatenate((X_plot, Xn), axis=1)

X_plot = np.matrix(X_plot).reshape(m, plot_degree * n + 1)
Y_plot = np.matrix(data_matrix[:, -1])


steepest_descent_deg_2, logs_deg_2 = steepest_descent(X_plot, Y_plot, initial_theta_plot, epochs = 60)
momentum_deg_2, logs_momentum_2 = momentum(X_plot, Y_plot, initial_theta_plot, epochs = 60, degree=1)

cost_history = [row[0] for row in logs_deg_2]
all_thetas = [row[1] for row in logs_deg_2]


theta0_history = [row[0].item() for row in all_thetas]
theta1_history = [row[1].item() for row in all_thetas]

cost_history = np.array(cost_history)
theta0_history = np.array(theta0_history)
theta1_history = np.array(theta1_history)

cost_history_momentum = [row[0] for row in logs_momentum_2]
all_thetas_momentum = [row[1] for row in logs_momentum_2]

theta0_history_momentum = [row[0].item() for row in all_thetas_momentum]
theta1_history_momentum = [row[1].item() for row in all_thetas_momentum]

cost_history_momentum = np.array(cost_history_momentum)
theta0_history_momentum = np.array(theta0_history_momentum)
theta1_history_momentum = np.array(theta1_history_momentum)

theta0_vals = theta0_history
theta1_vals = theta1_history
J_vals = np.zeros((len(theta0_vals), len(theta1_vals)))

c1=0
c2=0
pom = 0
for i in theta0_vals:
    for j in theta1_vals:
        t = np.array([i, j])
        J_vals[c1][c2] = cost_history[pom]
        c2=c2+1
    c1=c1+1
    pom += 1
    c2=0

theta0_vals_momentum = theta0_history_momentum
theta1_vals_momentum = theta1_history_momentum
J_vals = np.zeros((len(theta0_vals_momentum), len(theta1_vals_momentum)))

c1=0
c2=0
pom = 0
for i in theta0_vals_momentum:
    for j in theta1_vals_momentum:
        t = np.array([i, j])
        J_vals[c1][c2] = cost_history_momentum[pom]
        c2=c2+1
    c1=c1+1
    pom += 1
    c2=0

import plotly.express as px
fig = px.line_3d(x=theta0_vals, y=theta1_vals, z=[x[0] for x in J_vals])
fig.show()
fig2 = px.line_3d(x=theta0_vals_momentum, y=theta1_vals_momentum, z=[x[0] for x in J_vals])
fig2.show()

import plotly.graph_objects as go

fig = go.Figure(data=[go.Surface(x=theta0_vals, y=theta1_vals, z=J_vals)])
fig.update_layout(title='Loss function for different thetas', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

import plotly.graph_objects as go

fig = go.Figure(data=[go.Surface(x=theta0_vals_momentum, y=theta1_vals_momentum, z=J_vals)])
fig.update_layout(title='Loss function for different thetas', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

6. Zjawisko over- i underfittingu

Jednymi z najważniejszych podstawowych pojęć nauki o danych są pojęcia nadmiernego dopasowania i generalizacji. Jeśli pozwolimy sobie na wystarczającą elastyczność podczas poszukiwania wzorców w konkretnym zbiorze danych, to te wzorce znajdziemy. Niestety, te „wzorce” mogą być tylko zdarzeniami losowymi w ramach danych. W ogólności, chcemy aby model poprawnie przewidywał nowe dane, więc jest to zjawisko nieporządane.

Innym zjawiskiem jest niedostateczne dopasowanie. Jest to sytuacja, gdy model nie potrafi właściwie ani modelować danych treningowych, ani nie jest w stanie zrobić nic sensownego z nowymi danymi.

6. Regresja wielomianowa z wykorzystaniem gotowych bibliotek

from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import Ridge, LinearRegression

model = make_pipeline(PolynomialFeatures(degree=degree, include_bias=True), 
                       LinearRegression())
model.fit(data[["Height"]],Y)

Pipeline(steps=[('polynomialfeatures', PolynomialFeatures(degree=4)),
                ('linearregression', LinearRegression())])

Y_plot_2 = model.predict([[x] for x in X_plot_lib])

fig = plt.figure(figsize=(10,5))
chart = fig.add_subplot()
chart.plot(data["Height"], Y ,"go")
chart.plot(X_plot_lib, Y_plot_2, color="red", lw=2, label=f"degree {degree}")
plt.ylim([0,250])
plt.xlim([0,100])
degree

4

polyfit_theta = np.polyfit(data["Height"]/data["Height"].max(), Y, degree)
plot_polynomial_regression(polyfit_theta)