t2_circ.py

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

rng = np.random.RandomState(0)
M, n_sensors = 100, 1000
# Make "Fourier coefficients" here
data = rng.randn(M, n_sensors) + rng.randn(M, n_sensors) * 1j
data += 0.  # can be non-zero to test that it actually works for some signal
mean = np.mean(data, axis=0)
"""
fig, ax = plt.subplots(1, constrained_layout=True)
x = y = np.zeros(data.size)
ax.quiver(x, y, data.real.ravel(), data.imag.ravel(),
          scale=1, units='xy', alpha=0.1)
ax.quiver(0, 0, mean.real.mean(), mean.imag.mean(),
          scale=1, units='xy', color='r', lw=3)
ax.set(xlim=[-3, 3], ylim=[-3, 3])
"""
# eq3:
t2_circ = (M - 1) * np.abs(mean) ** 2 / np.sum(np.abs(data - mean) ** 2)
# M*T2circ is distributed according to F[2, 2M-2]
rv = stats.f(2, 2 * M - 2)
M_t2_circ = M * t2_circ
fig, ax = plt.subplots(1, figsize=(5, 3), constrained_layout=True)
bins = int(round(n_sensors / 30))
hist, bin_edges = np.histogram(M_t2_circ, bins=bins)
w = (bin_edges[1:] - bin_edges[:-1])
hist = hist / (hist * w).sum() / n_sensors
ax.step(bin_edges, np.concatenate([hist, [0]]), where='post')
ps = rv.pdf(bin_edges)
ax.plot(bin_edges, ps, color='orange')
# eq4:
sig = M_t2_circ > rv.isf(0.05)
print(sig.mean())  # should be ~0.05