Generate a data set by sampling 100 points from $y = \text{tanh}(x)$ between $x = -5$ to $5$, adding random noise with a magnitude of $0.25$ to $y$. Fit it with cluster-weighted modeling, using a local linear model for the clusters, $y = a + bx$. Plot the data and the resulting forecast $\langle y | x \rangle$, uncertainty $\langle \sigma_y^2 | x \rangle$, and cluster probabilities $p(\vec{x}|c_m)p(c_m)$ as the number of clusters is increased (starting with one), and animate their evolution during the EM iterations.
# From Neil's cwm.py code
import matplotlib.pyplot as plt
import numpy as np
def run(n_clusters):
x_min = -5
x_max = 5
n_pts = 100
min_var = 0.001
n_iterations = 100
x = np.arange(x_min, x_max, (x_max - x_min)/n_pts)
n_plot = n_pts
y = np.tanh(x) + 0.25*np.random.rand(n_pts)
x_plot = np.copy(x)
mux = x_min + (x_max - x_min)*np.random.rand(n_clusters)
beta = np.random.rand(2, n_clusters)
var_x = np.ones(n_clusters)
var_y = np.ones(n_clusters)
pc = np.ones(n_clusters)/n_clusters
fig,ax = plt.subplots()
plt.ion()
for step in range(n_iterations):
pxgc = np.exp(-(np.outer(x,np.ones(n_clusters))\
-np.outer(np.ones(n_pts), mux))**2/(2*np.outer(np.ones(n_pts), var_x))) \
/np.sqrt(2*np.pi*np.outer(np.ones(n_pts), var_x))
pygxc = np.exp(-(np.outer(y, np.ones(n_clusters)) \
- (np.outer(np.ones(n_pts), beta[0,:]) + np.outer(x, beta[1:])))**2 \
/(2*np.outer(np.ones(n_pts), var_y))) \
/np.sqrt(2*np.pi*np.outer(np.ones(n_pts), var_y))
pxc = pxgc*np.outer(np.ones(n_pts), pc)
pxyc = pygxc*pxc
pcgxy = pxyc/np.outer(np.sum(pxyc, 1), np.ones(n_clusters))
pc = np.sum(pcgxy, 0)/n_pts
mux = np.sum(np.outer(x, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc)
var_x = min_var + np.sum((np.outer(x, np.ones(n_clusters)) \
- np.outer(np.ones(n_pts), mux))**2*pcgxy, 0)/(n_pts*pc)
a = np.array([\
np.sum(np.outer(y*1, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc),\
np.sum(np.outer(y*x, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc)\
])
B = np.array([[\
np.ones(n_clusters),\
np.sum(np.outer(1*x, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc),\
],[\
np.sum(np.outer(x*1, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc),\
np.sum(np.outer(x*x, np.ones(n_clusters))*pcgxy, 0)/(n_pts*pc),\
]])
for c in range(n_clusters):
beta[:,c] = np.linalg.inv(B[:,:,c])@a[:,c]
var_y = min_var + np.sum((np.outer(y, np.ones(n_clusters))\
- (np.outer(np.ones(n_pts), beta[0,:]) + np.outer(x, beta[1:])))**2\
*pcgxy, 0)/(n_pts*pc)
plt.cla()
plt.plot(x, y, '.')
y_plot = np.sum((np.outer(np.ones(n_plot), beta[0,:])\
+ np.outer(x, beta[1,:]))*pxc, 1)/np.sum(pxc, 1)
plt.plot(x_plot, y_plot)
y_var_plot = np.sum((var_y + (np.outer(np.ones(n_plot), beta[0,:])\
+ np.outer(x, beta[1,:]))**2)*pxc, 1)/np.sum(pxc, 1) - y_plot**2
y_std_plot = np.sqrt(y_var_plot)
plt.plot(x_plot, y_plot + y_std_plot, 'k--')
plt.plot(x_plot, y_plot - y_std_plot, 'k--')
y_plot = 0.8*pxc/np.max(pxc) - 2
ax.plot(x_plot, y_plot, 'k')
fig.canvas.draw()
plt.show()
run(1)
run(2)
run(3)
run(4)
run(5)
run(6)