Origins of the Von Mises distribution
The goal here is to check if the Von Mises distribution is the a priori choice to make when handling polar coordinates.
Let's first initialize the notebook:
import numpy as np
np.set_printoptions(precision=6, suppress=True)
import os
%matplotlib inline
%config InlineBackend.figure_format='retina'
#%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
phi = (np.sqrt(5)+1)/2
fig_width = 8
figsize = (fig_width, fig_width/phi)
N = 100
m_x, m_y = 5., 3
sigma_x, sigma_y = 1., .2
x, y = m_x + sigma_x*np.random.randn(N), m_y + sigma_y*np.random.randn(N)
x.min(), x.max()
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width))
ax.scatter(x, y)
ax.set_xlim(0, 10)
ax.set_ylim(0, 10);
Let's do the same visualization but now transformed into polar coordinates:
r = np.sqrt(x**2 + y**2)
theta = np.arctan2(y, x)
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width), subplot_kw=dict(polar=True))
ax.scatter(theta % np.pi, r);
ax.set_ylim(0, 10);
projection on the ring: the von Mises distribution¶
Let's now do the histogram of these values in the polar space of angles:
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width))
ax.hist(theta, bins=np.linspace(0, np.pi, 10));
ax.set_xlabel('theta')
ax.set_ylabel('#');
#ax.set_ylim(0, 10);
Can we fit it with a Von Mises distribution?
N = 1000000
N_bins = 512
m_x, m_y = 5., 3
sigma_x, sigma_y = 1., .2
x, y = m_x + sigma_x*np.random.randn(N), m_y + sigma_y*np.random.randn(N)
r = np.sqrt(x**2 + y**2)
theta = np.arctan2(y, x)
bins = np.linspace(0, np.pi, N_bins)
v_hist, v_edges = np.histogram(theta, bins=bins, density=True);
v_hist.shape, v_edges.shape
def VonMises(theta, amp, theta0, Btheta):
return amp * np.exp((np.cos(2*(theta-theta0))-1) / 4 / Btheta**2)
from lmfit import Model
theta = .5*( bins[1:] + bins[:-1] )
vm_model = Model(VonMises)
result = vm_model.fit(v_hist, theta=theta, amp=1., theta0=1.5, Btheta=.5)
print(result.fit_report())
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width))
ax.plot(theta, v_hist, 'bo')
ax.plot(theta, result.init_fit, 'k--')
ax.plot(theta, result.best_fit, 'r-')
ax.set_xlabel('theta')
ax.set_ylabel('#');
Seems pretty good... Let's test for a larger set of Gaussian distributions with random parameters:
N_try = 14
fig, axs = plt.subplots(N_try, 1, figsize=(fig_width, fig_width))
np.random.seed(42)
redchis, Bthetas = [], []
M, S = .25, 1.
for ax in axs:
m_x, m_y = M*np.random.randn(), M*np.random.randn()
sigma_x, sigma_y = S*np.random.randn(), S*np.random.randn()
x, y = m_x + sigma_x*np.random.randn(N), m_y + sigma_y*np.random.randn(N)
thetas = np.arctan2(y, x)
v_hist, v_edges = np.histogram(thetas, bins=bins, density=True);
vm_model = Model(VonMises)
result = vm_model.fit(v_hist, theta=theta, amp=3., theta0=1.5, Btheta=.5)
#print(result.fit_report())
redchis.append(result.redchi)
Bthetas.append(np.abs(result.best_values['Btheta']))
ax.plot(theta, v_hist, 'o', ms=1)
ax.plot(theta, result.best_fit, '-')
ax.set_xlabel('theta')
ax.set_ylim(0, 4);
ax.set_ylabel('#');
print('redchis=', np.mean(redchis), ' +/-', np.std(redchis))
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width))
ax.scatter(Bthetas, redchis)
ax.set_ylim(0)
ax.set_xlabel('Btheta')
ax.set_ylabel('redchi');
Pretty good, but not perfect (especially for peaked distributions) ... If the distribution is centered ($m_x = m_y = 0$) we get similar (yet better) results:
N_try = 14
fig, axs = plt.subplots(N_try, 1, figsize=(fig_width, fig_width))
np.random.seed(42)
redchis, Bthetas = [], []
M, S = .0, 1.
for ax in axs:
m_x, m_y = M*np.random.randn(), M*np.random.randn()
sigma_x, sigma_y = S*np.random.randn(), S*np.random.randn()
x, y = m_x + sigma_x*np.random.randn(N), m_y + sigma_y*np.random.randn(N)
thetas = np.arctan2(y, x)
v_hist, v_edges = np.histogram(thetas, bins=bins, density=True);
vm_model = Model(VonMises)
result = vm_model.fit(v_hist, theta=theta, amp=3., theta0=1.5, Btheta=.5)
#print(result.fit_report())
redchis.append(result.redchi)
Bthetas.append(np.abs(result.best_values['Btheta']))
ax.plot(theta, v_hist, 'o', ms=1)
ax.plot(theta, result.best_fit, '-')
ax.set_xlabel('theta')
ax.set_ylim(0, 4);
ax.set_ylabel('#');
print('redchis=', np.mean(redchis), ' +/-', np.std(redchis))
fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width))
ax.scatter(Bthetas, redchis)
ax.set_ylim(0)
ax.set_xlabel('Btheta')
ax.set_ylabel('redchi');
Yet, in Wang, F., and Gelfand, A. E. 2013. Directional data analysis under the general projected normal distribution. Statistical methodology 10(1):113–127., it is claimed that :
The von Mises distribution is the maximum entropy distribution on the circle for a fixed mean direction and dispersion and can be considered as an analogue of the normal distribution on the real line.
Still, I have nbo proof for that - this would certainly begin with a change of variables
Is there something better?
some book keeping for the notebook¶
%load_ext watermark
%watermark -i -h -m -v -p numpy,matplotlib,lmfit -r -g -b