(16.1)

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.

In [10]:
# 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)
In [11]:
run(2)
In [12]:
run(3)
In [13]:
run(4)
In [14]:
run(5)
In [15]:
run(6)