In [1]:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
import matplotlib.colors as colors
import matplotlib as mpl
from matplotlib import rc
rc('animation', html='html5')
from matplotlib import collections
from matplotlib.animation import FuncAnimation
from ipywidgets import interact, interactive, fixed, interact_manual, IntSlider
In [2]:
pal = sns.color_palette()
def tand(x):
return np.tan(x * 3.141592653589793 / 180)
def sind(x):
return np.sin(x * 3.141592653589793 / 180)
def cosd(x):
return np.cos(x * 3.141592653589793 / 180)
In [3]:
np.random.seed(22)
mat_obs = np.random.multivariate_normal([0, 0], [[2, 0.7], [0.7, 1]], size=50)
mat_obs = mat_obs - mat_obs.mean(axis=0)
In [4]:
class InteractiveProjection(object):
def __init__(self, obs, model_type="PCA"):
self.obs = obs - obs.mean(axis=0)
self.model_type = model_type
self.vcv_obs = self.obs.T.dot(self.obs) # Sum of Square Matrix
eigval, eigvec = np.linalg.eigh(self.vcv_obs) # Eigenvalues and Eigenvectors
idx = eigval.argsort()[::-1]
self.eigval = eigval[idx]
self.eigvec = eigvec[idx,:]
self.angle = 0
self.minerror = 100000
self.vec_proj, self.obs_proj, self.error = self.update_projection(self.angle)
self.fig, self.ax, self.proj_dots, self.linecoll, self.proj_x, self.proj_y, self.annot1, self.annot2 = self.init_figure()
def init_figure(self):
eigvec, eigval = self.eigvec, self.eigval
obs = self.obs
pal = sns.color_palette()
a1 = -5.3
a2 = -4.3
fig, ax = plt.subplots(1, 1, figsize=(5, 5), sharey=True)
plt.autoscale(False)
sns.despine()
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.scatter(0, 0, s=50, color="white", zorder=100, edgecolors="black", linewidths=2)
ax.scatter(*obs.T, color=pal[1], s=20)
e1, = ax.plot(eigvec[0, 0]*np.array([-a2, -a1]), eigvec[0, 1]*np.array([-a2, -a1]), color=pal[4], lw=2)
e2, = ax.plot(eigvec[1, 0]*np.array([-a2, -a1]), eigvec[1, 1]*np.array([-a2, -a1]), color=pal[3], lw=2)
ax.plot(eigvec[0, 0]*np.array([a2, a1]), eigvec[0, 1]*np.array([a2, a1]), color=pal[4], lw=2)
ax.plot(eigvec[1, 0]*np.array([a2, a1]), eigvec[1, 1]*np.array([a2, a1]), color=pal[3], lw=2)
linecoll = collections.LineCollection([], color=pal[2], lw=1, linestyle=":", alpha=0.7)
ax.add_collection(linecoll)
annot1 = ax.annotate("\n", (1, -2.5), va="center")
annot2 = ax.annotate("\n", (-3.5, 2.5), va="center")
proj_dots = ax.scatter([], [], color=pal[2], s=10, animated=True)
proj_x, = ax.plot([], [], 'k', lw=1, alpha=0.7)
proj_y, = ax.plot([], [], 'k', lw=1, alpha=0.7)
plt.figlegend((e1, e2), ("First Eigenvector", "Second Eigenvector"), loc="lower center", ncol=2,
bbox_to_anchor=(0.5, 0), title="")
plt.tight_layout(rect=(0, 0.1, 1, 1))
plt.close()
return fig, ax, proj_dots, linecoll, proj_x, proj_y, annot1, annot2
def update_projection(self, angle):
obs = self.obs
# Angle of rotation
self.angle = angle
# Projections corresponding to angles
vec_proj = np.array([cosd(angle), sind(angle)]).reshape(2, -1)
# Projection matrices
if self.model_type == "PCA":
mat_proj = vec_proj.dot(vec_proj.T)
elif self.model_type == "OLSX":
mat_proj = np.array([[0, 0], [cosd(angle)/sind(angle), 1]])
else:
mat_proj = np.array([[1, sind(angle)/cosd(angle)], [0, 0]])
# Projected observations
obs_proj = obs.dot(mat_proj)
error = obs_proj.dot([[0, 1], [-1, 0]]) * (obs-obs_proj)
return vec_proj, obs_proj, error
def plot_projection(self, angle):
self.vec_proj, self.obs_proj, self.error = self.update_projection(angle)
obs = self.obs
w = self.vec_proj
obs_proj = self.obs_proj
self.proj_dots.set_offsets(obs_proj)
segs = np.rollaxis(np.dstack([obs, obs_proj]), 2, 1)
self.linecoll.set_segments(segs)
self.proj_x.set_data(w[0]*4.5*[-1, 1], w[1]*4.5*[-1, 1])
self.proj_y.set_data(-w[1]*1.75*[-1, 1], w[0]*1.75*[-1, 1])
error = np.square(self.error).sum()
self.annot1.set_text(f"Current Error:\n{error:.1f}")
if error < self.minerror:
self.minerror = error
self.bestangle = self.angle
self.annot2.set_text(f"Min Error So Far:\n{self.minerror:.1f} at {self.bestangle:.0f}°")
return self.fig
In [5]:
spring = InteractiveProjection(mat_obs, model_type="OLSX");
interact(spring.plot_projection, angle=IntSlider(min=-30,max=75,step=1,value=0));
In [6]:
class AnimationProjection(object):
def __init__(self, obs):
self.obs = obs - obs.mean(axis=0) # Mean-center, since rod is fixed at 0.
self.vcv_obs = self.obs.T.dot(self.obs) # Sum of Square Matrix
eigval, eigvec = np.linalg.eigh(self.vcv_obs) # Eigenvalues and Eigenvectors
idx = eigval.argsort()[::-1]
self.eigval = eigval[idx]
self.eigvec = eigvec[idx,:]
self.angles = {"PCA": -30,
"OLS1": -30,
"OLS2": -30}
self.vec_proj, self.obs_proj, self.error, self.energy = self.init_system()
self.fig, self.axes, self.proj_dots, self.linecoll, self.proj_x, self.proj_y, self.annot = self.init_figure()
return None
def init_figure(self):
eigvec, eigval = self.eigvec, self.eigval
obs = self.obs
pal = sns.color_palette()
a1 = -5.3
a2 = -4.3
titles = ["PCA\n(Min. Euclidian Distance)",
"X regressed on Y\n(Min. Horizontal Distance)",
"Y regressed on X\n(Min. Vertical Distance)"]
fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
plt.autoscale(False)
for ax, title in zip(axes, titles):
ax.set_title(title)
sns.despine()
ax.set_aspect("equal")
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.scatter(0, 0, s=50, color="white", zorder=100, edgecolors="black", linewidths=2)
ax.scatter(*obs.T, color=pal[1], s=20)
e1, = ax.plot(eigvec[0, 0]*np.array([-a2, -a1]), eigvec[0, 1]*np.array([-a2, -a1]), color=pal[4], lw=2)
e2, = ax.plot(eigvec[1, 0]*np.array([-a2, -a1]), eigvec[1, 1]*np.array([-a2, -a1]), color=pal[3], lw=2)
ax.plot(eigvec[0, 0]*np.array([a2, a1]), eigvec[0, 1]*np.array([a2, a1]), color=pal[4], lw=2)
ax.plot(eigvec[1, 0]*np.array([a2, a1]), eigvec[1, 1]*np.array([a2, a1]), color=pal[3], lw=2)
linecoll = {}
proj_dots = {}
proj_x = {}
proj_y = {}
annot = {}
for ax, label in zip(axes, ["PCA", "OLS1", "OLS2"]):
lines = collections.LineCollection([], color=pal[2], lw=1, linestyle=":", alpha=0.7)
linecoll[label] = lines
ax.add_collection(lines)
annot[label] = ax.annotate("", (0, -3))
proj_dots[label] = ax.scatter([], [], color=pal[2], s=10, animated=True)
proj_x[label], = ax.plot([], [], 'k', lw=1, alpha=0.7)
proj_y[label], = ax.plot([], [], 'k', lw=1, alpha=0.7)
plt.figlegend((e1, e2), ("First Eigenvector", "Second Eigenvector"), loc="lower center", ncol=2,
bbox_to_anchor=(0.5, 0), title="")
plt.tight_layout(rect=(0, 0.1, 1, 1))
return fig, axes, proj_dots, linecoll, proj_x, proj_y, annot
def init_system(self):
obs = self.obs
# Angle of rotation
angles = self.angles
# Projections corresponding to angles
vec_proj = {k:np.array([cosd(v), sind(v)]).reshape(2, -1) for k, v in angles.items()}
# Projection matrices
mat_proj = {
"PCA": vec_proj["PCA"].dot(vec_proj["PCA"].T),
"OLS1": np.array([[0, 0], [cosd(angles["OLS1"])/sind(angles["OLS1"]), 1]]),
"OLS2": np.array([[1, sind(angles["OLS2"])/cosd(angles["OLS2"])], [0, 0]])
}
# Projected observations
obs_proj = {k:obs.dot(v) for k, v in mat_proj.items()}
error = {k:v.dot([[0, 1], [-1, 0]]) * (obs-v) for k, v in obs_proj.items()}
energy = {k: v.sum() for k, v in error.items()}
return vec_proj, obs_proj, error, energy
def update_system(self):
energy = {k: (v + self.error[k].sum())*0.93 for k, v in self.energy.items()}
angles = {k: (v + self.energy[k]/40) for k, v in self.angles.items()}
obs = self.obs
# Projections corresponding to angles
vec_proj = {k:np.array([cosd(v), sind(v)]).reshape(2, -1) for k, v in angles.items()}
# Projection matrices
mat_proj = {
"PCA": vec_proj["PCA"].dot(vec_proj["PCA"].T),
"OLS1": np.array([[0, 0], [cosd(angles["OLS1"])/sind(angles["OLS1"]), 1]]),
"OLS2": np.array([[1, sind(angles["OLS2"])/cosd(angles["OLS2"])], [0, 0]])
}
# Projected observations
obs_proj = {k:obs.dot(v) for k, v in mat_proj.items()}
error = {k:v.dot([[0, 1], [-1, 0]]) * (obs-v) for k, v in obs_proj.items()}
self.angles, self.vec_proj, self.obs_proj, self.error, self.energy = angles, vec_proj, obs_proj, error, energy
return None
def init_projection(self):
for ax, label in zip(self.axes, ["PCA", "OLS1", "OLS2"]):
obs = self.obs
w = self.vec_proj[label]
obs_proj = self.obs_proj[label]
self.proj_dots[label].set_offsets(obs_proj)
segs = np.rollaxis(np.dstack([obs, obs_proj]), 2, 1)
self.linecoll[label].set_segments(segs)
self.proj_x[label].set_data(w[0]*4.5*[-1, 1], w[1]*4.5*[-1, 1])
if label == "PCA":
self.proj_y[label].set_data(-w[1]*1.75*[-1, 1], w[0]*1.75*[-1, 1])
self.annot[label].set_text("Current Error:\n{:.1f}".format(np.square(self.error[label]).sum()))
return self.linecoll[label],
def plot_projection(self, i):
self.update_system()
for ax, label in zip(self.axes, ["PCA", "OLS1", "OLS2"]):
obs = self.obs
w = self.vec_proj[label]
obs_proj = self.obs_proj[label]
self.proj_dots[label].set_offsets(obs_proj)
segs = np.rollaxis(np.dstack([obs, obs_proj]), 2, 1)
self.linecoll[label].set_segments(segs)
self.proj_x[label].set_data(w[0]*4.5*[-1, 1], w[1]*4.5*[-1, 1])
if label == "PCA":
self.proj_y[label].set_data(-w[1]*1.75*[-1, 1], w[0]*1.75*[-1, 1])
self.annot[label].set_text("Current Error:\n{:.1f}".format(np.square(self.error[label]).sum()))
return self.linecoll[label],
def animate(self, frames, interval):
anim = FuncAnimation(self.fig, init_func=self.init_projection,
func=self.plot_projection, frames=frames, blit=True, interval=interval)
return anim
In [7]:
spring = AnimationProjection(mat_obs)
plt.close()
spring.animate(500, 30)
Out[7]:
In [ ]: