In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
np.random.seed(17)
In [3]:
def poly(a):
return lambda x: np.dot(a, [x ** i for i in range(len(a))])
In [4]:
def apply(f, X):
return np.array([f(x) for x in X])
In [5]:
def sub_plot(ax, x, y, n, bounds, *methods, **params):
gs = [method(x, y, n, **params) for method in methods]
points = np.linspace(*bounds, num=100)
for i in range(len(methods)):
ax.plot(points, apply(gs[i], points), label=methods[i].__name__)
ax.plot(x, y, 'rx')
title = " vs ".join([method.__name__.capitalize() for method in methods]) + " approximation. " + \
str(n) + " degree."
if len(params) > 0:
title += " Params: " + str(params)
ax.set_title(title)
ax.legend(loc="best")
In [6]:
def plot(x, y, n, bounds, *methods, **params):
figure, ax = plt.subplots(figsize=(15, 10))
sub_plot(ax, x, y, n, bounds, *methods, **params)
In [7]:
def ols(x, y, n, verbose=False, **unused_params):
H = np.zeros((len(x), n + 1))
for i in range(len(x)):
for j in range(n + 1):
H[i][j] = x[i] ** j
A = np.dot(H.T, H)
z = np.dot(H.T, y)
c = np.linalg.solve(A, z)
if not verbose:
print(c)
return poly(c)
In [8]:
x = np.linspace(-np.pi / 2.1, np.pi / 2.1, num=50)
y = np.tan(x) + np.random.normal(0, 0.25, 50)
In [9]:
plot(x, y, 3, (-2, 2), ols)
In [22]:
x = np.linspace(0, np.pi, num=100)
y = np.sin(x**2) + np.random.normal(0, 0.05, 100)
In [23]:
plot(x, y, 4, (-0.5, 3.5), ols)
In [12]:
def n_plot(x, y, ns, bounds, *methods):
figure, axs = plt.subplots(len(ns), figsize=(15, len(ns) * 5))
for i in range(len(ns)):
sub_plot(axs[i], x, y, ns[i], bounds, *methods)
In [13]:
x = np.linspace(-1, 3 * np.pi, num=50)
y = np.sin(x)*x + np.random.normal(0, 0.2, 50)
In [14]:
n_plot(x, y, np.arange(1, 5), (-2, 10), ols)
In [15]:
def smoothing_spline(x, y, n, alpha=0.1, verbose=False):
H = np.zeros((len(x), n + 1))
for i in range(len(x)):
for j in range(n + 1):
H[i][j] = x[i] ** j
D = np.zeros((n + 1, n + 1))
for i in range(n + 1):
for j in range(n + 1):
if j == i or j == i + 2:
D[i][j] = 1
elif j == i + 1:
D[i][j] = -2
A = np.dot(H.T, H) + alpha * np.dot(D.T, D)
z = np.dot(H.T, y)
c = np.linalg.solve(A, z)
if not verbose:
print(c)
return poly(c)
In [16]:
x = np.linspace(0, np.pi, num=100)
y = np.arctan(x**8) + np.random.normal(0, 0.05, 100)
In [17]:
plot(x, y, 5, (-0.5, 3.5), smoothing_spline, ols, alpha=0.1)
In [18]:
def n_plot(x, y, n, alphas, bounds, *methods):
figure, axs = plt.subplots(len(alphas), figsize=(15, len(alphas) * 5))
for i in range(len(alphas)):
sub_plot(axs[i], x, y, n, bounds, *methods, alpha=alphas[i])
In [19]:
x = np.linspace(0.01, np.pi, num=20)
y = np.cos(x**3)*x**0.1 + np.random.normal(0, 0.1, 20)
In [20]:
n_plot(x, y, 9, np.logspace(-1, 2, num=4, endpoint=True), (-0.5, 3.2), smoothing_spline)