nonlinear_gd.py

"""
@classAndrew
Nonlinear regression with a single layer perceptron.
Forgot when this was written
"""

import numpy as np
import statistics
import matplotlib.pyplot as plt

PREC = 500
LEFT = -10
RIGHT = 10
fuzz = (np.random.rand(PREC) - 0.5)*0.8

N = 7
#significance = [10**-i for i in range(N-1, -1, -1)] # this is for polynomials
rgen = np.random.RandomState(27)
coef = (rgen.rand(N)-0.5)*5#*significance

x = np.arange(LEFT, RIGHT, (RIGHT-LEFT)/PREC)

def F_n(x, coef):
    return coef[0]+sum(coef[i]*np.sin(coef[i+1]*x) for i in range(1, len(coef), 2))

class StandardScaler:
    def __init__(self):
        self.mean = 0
        self.stdev = 0
    
    def fit_transform(self, X):
        self.mean = statistics.mean(X)
        self.stdev = statistics.stdev(X)
        return (X-self.mean)/self.stdev
    
    def inverse_transform(self, X):
        return X*self.stdev + self.mean

class NonlinearGD:
    def __init__(self, N, P_n, w_scale = 1, eta=0.001, max_iter=2000):
        self.eta = eta
        self.max_iter = max_iter
        self.N = N
        self.w_ = None
        self.P_n = P_n
        self.w_scale = w_scale
        self.cost_ = []

    def fit(self, X, y):
        self.w_ = rgen.normal(loc=0., scale=self.w_scale, size=1 + self.N)

        for _ in range(self.max_iter):
            # print(dJ_dw)
            self.w_ -= self.eta*self.grad(X, y)

        return self

    def predict(self, X): # -> vector
        # transform X to a len(X) by N (exclude bias) matrix
        return self.P_n(X, self.w_[1:])
    
    def net_input(self, X):
        return self.P_n(X, self.w_[1:]) # + self.w_[0]
    
    def sse(self, errors):
        return sum(errors**2)

    def grad(self, X, y):
        dw = 1e-5
        dJ_dw = np.zeros(len(self.w_))

        old = self.sse(y-self.net_input(X))
        self.cost_.append(old)

        for i in range(1, len(self.w_)):
            self.w_[i] += dw
            dJ_dw[i] = (self.sse(y-self.net_input(X))-old)/dw
            self.w_[i] -= dw

        return dJ_dw

y = F_n(x, coef)+fuzz
# print(pgd.w_scale)

def plot_scaled():
    pgd = NonlinearGD(N, F_n, w_scale=statistics.stdev(y)*2.5, eta=0.001989, max_iter=1500)
    y_sc = StandardScaler()
    x_sc = StandardScaler()

    y_std = y_sc.fit_transform(y)
    x_std = x_sc.fit_transform(x)
    pgd.fit(x_std, y_std)
    plt.plot(x, y_sc.inverse_transform(pgd.predict(x_std)), color="red", linewidth=2)
    print(f"Scaler SSE: {pgd.cost_[0]:.4f}, {pgd.cost_[-1]:.4f}")

def plot():
    pgd = NonlinearGD(N, F_n, w_scale=statistics.stdev(y)*2.5)
    pgd.fit(x, y)
    plt.plot(x, pgd.predict(x), color="lime", linewidth=2)
    print(f"Regular SSE: {pgd.cost_[0]:.4f}, {pgd.cost_[-1]:.4f}")

# plot()
plot_scaled()
# print(pgd.cost_[:10], "...", pgd.cost_[-10:])
# print(pgd.w_)

plt.scatter(x, y)
plt.show()