Demo: Motion illusions as optimal percepts

Weiss, Y., Simoncelli, E. P., & Adelson, E. H. (2002). Motion illusions as optimal percepts. Nature Neuroscience, 5(6), 598–604. doi:10.1038/nn858

First, just create our imports and define a few helper functions to get started:

%matplotlib inline
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from ipywidgets import interact

def imshow(ax, p):
    """Show the probabilities as a function of x and y velocities."""
    ax.imshow(p.T, origin='lower', interpolation='nearest', cmap='gray')
    xmid = (p.shape[1] - 1) / 2
    ymid = (p.shape[0] - 1) / 2
    ax.vlines(ymid, 0, p.shape[1], color='gray')
    ax.hlines(xmid, 0, p.shape[0], color='gray')
    ax.set_xticks([])
    ax.set_yticks([])

def uniform(x, low, high):
    """Compute the log probability for a uniform random variable between (low, high)."""
    return scipy.stats.uniform.logpdf(x, low, high - low)

def norm(x, mu, sigma):
    """Compute the log probability for a Gaussian random variable with parameters μ and σ."""
    return scipy.stats.norm.logpdf(x, mu, sigma)

Define a few options for the prior. In the paper, they used the equivalent of prior1, but I’m also interested in comparing to a uniform prior and a Gaussian prior with different mean:

def prior1(vx, vy):
    """Zero-mean Gaussian prior with σ=25"""
    return norm(vx, 0, 25) + norm(vy, 0, 25)

def prior2(vx, vy):
    """Velocity Average (VA) Gaussian prior with σ=5"""
    return norm(vx, 17.88461538, 5) + norm(vy, -14.42307692, 5)

def prior3(vx, vy):
    """Uniform prior between -50 and 50"""
    return uniform(vx, -50, 100) + uniform(vy, -50, 100)

Now, define the full model. This assumes a thin rhombus, but the prior function and the contrast (i.e., inverse sigma) can be modified:

def model(prior_func, sigma):
    vx, vy = np.ogrid[-50:51,-50:51]
    prior = prior_func(vx, vy)
    lh1 = uniform(vx, -50, 100) + norm(vy - vx + 30, 0, sigma)
    lh2 = uniform(vx, -50, 100) + norm(vy - 1.5 * vx + 45, 0, sigma)
    posterior = lh1 + lh2 + prior
    MAP = np.unravel_index(np.argmax(posterior), posterior.shape)

    fig, axes = plt.subplots(1, 4)
    fig.set_size_inches(12, 3)

    imshow(axes[0], np.exp(prior))
    axes[0].set_title('Prior')
    
    imshow(axes[1], np.exp(lh1))
    axes[1].set_title('Likelihood 1')
    
    imshow(axes[2], np.exp(lh2))
    axes[2].set_title('Likelihood 2')
    
    imshow(axes[3], np.exp(posterior))
    axes[3].set_title('Posterior')
    axes[3].autoscale(enable=False)
    axes[3].plot(MAP[0], MAP[1], 'ro')