import matplotlib.pyplot as plt
import numpy as np

from neurodsp.sim import sim_oscillation
from timescales.sim import sim_ar, sim_ar_spectrum

from timescales.conversions import phi_to_tau, tau_to_knee
from timescales.fit import PSD
from timescales.plts import set_default_rc
set_default_rc()

05. Estimation Bias

Bias timescale estimates arise from oscillatory activity and finite bias.

Oscilllations

To account for oscillatory bias, it is recommended to used the FOOOF model to explicitly account for spectral peaks or robust regression loss functions, such as huber or cauchy. Below we show how spectral peaks influence aperiodic fits when spectral peaks are not accounted for.

# Settings
n_seconds = 100
fs = 1000
f_range = (1, 100)
phi = 0.95
osc_freq = 10

# Simulate
sig_ap = sim_ar(n_seconds, fs, phi)
sig_pe = sim_oscillation(n_seconds, fs, osc_freq, mean=0, variance=.5)
sig = sig_ap + sig_pe

# Fit: robust regression
psd_robust = PSD()
psd_robust.compute_spectrum(sig, fs, f_range=f_range)
psd_robust.fit(method='cauchy')

psd_linear = PSD()
psd_linear.compute_spectrum(sig, fs, f_range=f_range)
psd_linear.fit(method='linear')

plt.figure(figsize=(10, 8))

plt.loglog(psd_robust.freqs, psd_robust.powers, label='Power', alpha=.5, lw=3)
plt.loglog(psd_robust.freqs, psd_robust.powers_fit, label='Robust Fit', lw=3)
plt.loglog(psd_linear.freqs, psd_linear.powers_fit, label='Linear Fit', lw=3)

plt.title('Linear vs Robust Regression')
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.legend();

Finite

Finite bias may be account for using autoregressive power spectral density instead of non-parametric estimates, such as Welch’s method. Accounting for finite bias allows the use of short time windows to explore the time or task-dependent nature of timescales.

# Settings
n_seconds = 2
fs = 1000

# Simulate
np.random.seed(0)
sig = sim_ar(n_seconds, fs, phi)

# AR
psd_ar = PSD()
psd_ar.compute_spectrum(sig, fs, ar_order=100)
psd_ar.fit(method="ar", f_scale=1.0)

plt.figure(figsize=(10, 8))

plt.loglog(psd_ar.freqs[1:], psd_ar.powers[1:],
           label='AR(100)', lw=3)

plt.loglog(psd_ar.freqs[1:], psd_ar.powers_fit[1:],
           label='AR(1) Fit\nof AR(100)', lw=3)

powers_true = sim_ar_spectrum(psd_ar.freqs, fs, phi, 1/fs)
plt.loglog(psd_ar.freqs[1:], powers_true[1:],
           label='Truth', lw=3, ls='--', color='k')

plt.title('Finite Bias')
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.legend();

Oscillations + Finite

Autoregressive PSD estimates preserve oscillatory structure by increasing the order of the model.

# Settings
n_seconds = 1
osc_freq = 10
ar_order = 50

# Simulate
np.random.seed(0)
sig = sim_ar(n_seconds, fs, phi)
sig += sim_oscillation(n_seconds, fs, osc_freq, variance=1)

# Fit
psd_ar = PSD()
psd_ar.compute_spectrum(sig, fs, ar_order=ar_order, f_range=(1e-3, np.inf))
psd_ar.fit(method='ar', f_scale=.2)

# Plot
plt.figure(figsize=(10, 8))

plt.loglog(psd_ar.freqs, psd_ar.powers,
           lw=3, color='C0', label="AR(100)")

plt.loglog(psd_ar.freqs, psd_ar.powers_fit,
           lw=3, color='C1', label="AR(1) Fit")

powers_true = sim_ar_spectrum(psd_ar.freqs, fs, phi, 1/fs)
plt.loglog(psd_ar.freqs, powers_true, color='k', ls='--', label='Truth')

plt.title('Finite + Oscillatory Bias')
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.legend();