NMM 2020 Jake Read 10 Function Architectures: Fitting¶
In [1]:
import numpy as np
import math
import matplotlib.pyplot as plt
import scipy as sp
14.2¶
Take as a data set $x = \{-10, -9, ..., 9, 10\}$ and $y(x) = 0$ if $x \leq 0$ and $y(x) = 1$ if $x > 0$
In [2]:
def y(x):
if(x <= 0):
return 0
else:
return 1
t = np.arange(-10, 11, 1)
d = np.zeros(len(t))
for i in range(len(t)):
d[i] = y(t[i])
plt.plot(t,d)
Out[2]:
a) fit the data with polynomials of 5, 10, and 15 terms, using a psuedo-inverse of the Vandermonde matrix.
In [3]:
def poly(t, p): # evaluate polynomial at t, with parameters p
sum = 0
for i in range(len(p)):
sum += p[i]*t**i
return sum
def polyseries(ts, p): # call the above for an array of t's
res = np.zeros(len(ts))
for i in range(len(ts)):
res[i] = poly(ts[i], p)
return res
def polyfactors(t, p): # return coefficients for the parameters p, at time t
res = np.zeros(len(p))
for i in range(len(p)):
res[i] = t**i
return res
t0 = np.arange(-10, 10, 0.001)
d0 = polyseries(t0, np.ones(5))
plt.plot(t0, d0)
Out[3]:
ok, now we need to find the p that best approximates the step function, with n+1 points we can exactly fit a polynomial with n components, so any less (5) will under-fit and (15) will over-fit, is my guess
In [4]:
#need to construct the vandermonde matrix
r0 = polyfactors(5, np.zeros(5))
print(r0)
#one of those, for each of our samples: first for 5, then 10, 15
vm5 = np.zeros((len(t), 5))
for i in range(len(t)):
#print(t[i])
vm5[i] = polyfactors(t[i], np.zeros(5))
vm5 # matrix for 5 terms,
Out[4]:
In [5]:
px5, r, rank, s = np.linalg.lstsq(vm5, d, rcond=None)
px5
Out[5]:
In [6]:
pd5 = polyseries(t0, px5)
plt.plot(t0, pd5)
plt.plot(t, d)
Out[6]:
doesn't seem radically unlikely that that's a best fit with just 5 terms, I'll try 10 and 15, to see if they improve.
In [7]:
vm10 = np.zeros((len(t), 10))
for i in range(len(t)):
vm10[i] = polyfactors(t[i], np.zeros(10))
# for reasons such as the above I will make no claim to code elegance
px10, r, rank, s = np.linalg.lstsq(vm10, d, rcond=None)
In [34]:
pd10 = polyseries(t0, px10)
plt.plot(t0, pd10)
for i in range(21):
plt.axvline(i-10) # visualize on-sample and off errors,
plt.plot(t, d)
Out[34]:
In [9]:
vm15 = np.zeros((len(t), 15))
for i in range(len(t)):
vm15[i] = polyfactors(t[i], np.zeros(15))
px15, r, rank, s = np.linalg.lstsq(vm15, d, rcond=None)
pd15 = polyseries(t0, px15)
plt.plot(t0, pd15)
plt.plot(t, d)
Out[9]:
so, that 15 term is a mystery to me. is this meant to be a discovery? the x^15 term is hard to fight off, if I manually clip it...
In [31]:
pd15 = polyseries(t0, px15[0:10])
plt.plot(t0, pd15)
plt.plot(t, d)
Out[31]:
b) do the same, but using RBFs of $r^3$
In [11]:
def r3(t, offset):
return (t+offset)**3
def r3factors(t, offsets):
res = np.zeros(len(offsets))
for i in range(len(offsets)):
res[i] = r3(t, offsets[i])
return res
def r3seq(offsets, oa, t): # offsets are positions, oa are coefficients for this offset
sum = 0
for i in range(len(offsets)):
sum += r3(t, offsets[i])*oa[i]
return sum
offsets5 = np.arange(-10, 10.1, 5)
oa5 = np.ones(5)
of5ex = np.zeros(len(t0))
for i in range(len(t0)):
of5ex[i] = r3seq(offsets5, oa5, t0[i])
plt.plot(t0, of5ex)
Out[11]:
In [12]:
# ahn matrix for 5,
of5 = np.zeros((len(t), 5))
for i in range(len(t)):
#print(t[i])
of5[i] = r3factors(t[i], offsets5)
of5
Out[12]:
In [13]:
ox5, r, rank, s = np.linalg.lstsq(of5, d, rcond=None)
ox5
Out[13]:
In [14]:
# reconstruct,
oxd5 = np.zeros(len(t0))
for i in range(len(t0)):
oxd5[i] = r3seq(offsets5, ox5, t0[i])
plt.plot(t0, oxd5)
plt.plot(t, d)
Out[14]:
In [15]:
# ok,
offsets10 = np.linspace(-10, 10, 10)
of10 = np.zeros((len(t), 10))
for i in range(len(t)):
#print(t[i])
of10[i] = r3factors(t[i], offsets10)
ox10, r, rank, s = np.linalg.lstsq(of10, d, rcond=None)
oxd10 = np.zeros(len(t0))
for i in range(len(t0)):
oxd10[i] = r3seq(offsets10, ox10, t0[i])
plt.plot(t0, oxd10)
plt.plot(t, d)
Out[15]:
In [16]:
# ok,
offsets15 = np.linspace(-10, 10, 15)
of15 = np.zeros((len(t), 15))
for i in range(len(t)):
#print(t[i])
of15[i] = r3factors(t[i], offsets15)
ox15, r, rank, s = np.linalg.lstsq(of15, d, rcond=None)
oxd15 = np.zeros(len(t0))
for i in range(len(t0)):
oxd15[i] = r3seq(offsets15, ox15, t0[i])
plt.plot(t0, oxd15)
plt.plot(t, d)
Out[16]:
In [17]:
ox15
Out[17]:
In [18]:
ox10
Out[18]:
In [19]:
ox5
Out[19]:
it's unclear to me why these all appear identical, but it feels like composing a step function with lots of little step functions is not an awesome idea, so maybe this is meant to be broken. let's see about off-axis errors
In [21]:
ts = np.linspace(-10.5, 10.5, 22)
ts
Out[21]:
In [ ]: