Transfer-Function and Prediction Diagnostics¶

These diagnostics answer different questions about one fitted model. Transfer-function plots describe the learned mapping; coherence and cross-spectra compare predictions with observations and should therefore be computed on held-out data.

In [1]:

import matplotlib
matplotlib.use("module://matplotlib_inline.backend_inline")
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import fftconvolve

from fftrf import TRF, pearsonr, r2_score

import matplotlib matplotlib.use("module://matplotlib_inline.backend_inline") import matplotlib.pyplot as plt import numpy as np from scipy.signal import fftconvolve from fftrf import TRF, pearsonr, r2_score

In [2]:

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12):
    rng = np.random.default_rng(seed)
    fs = 128.0
    tmin, tmax = 0.0, 0.25
    times = np.arange(round(tmin * fs), round(tmax * fs)) / fs
    true_kernel = (
        0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2)
        - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2)
    )

    stimulus, response = [], []
    for _ in range(n_trials):
        x = rng.standard_normal((n_samples, 1))
        y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples]
        y += noise_scale * rng.standard_normal(n_samples)
        stimulus.append(x)
        response.append(y[:, np.newaxis])
    return fs, tmin, tmax, times, true_kernel, stimulus, response

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12): rng = np.random.default_rng(seed) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmin * fs), round(tmax * fs)) / fs true_kernel = ( 0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2) - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2) ) stimulus, response = [], [] for _ in range(n_trials): x = rng.standard_normal((n_samples, 1)) y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples] y += noise_scale * rng.standard_normal(n_samples) stimulus.append(x) response.append(y[:, np.newaxis]) return fs, tmin, tmax, times, true_kernel, stimulus, response

Fit a model and reserve one trial for diagnostics¶

In [3]:

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials(
    seed=55,
    noise_scale=0.2,
)
train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

model = TRF(direction=1)
model.train(
    train_x, train_y,
    fs=fs, tmin=tmin, tmax=tmax,
    regularization=1e-2,
    segment_duration=1.0,
    overlap=0.5,
    window="hann",
)
_, heldout_r = model.predict(test_x, test_y)
print(f"held-out r: {float(heldout_r):.3f}")

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials( seed=55, noise_scale=0.2, ) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1] model = TRF(direction=1) model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_duration=1.0, overlap=0.5, window="hann", ) _, heldout_r = model.predict(test_x, test_y) print(f"held-out r: {float(heldout_r):.3f}")

held-out r: 0.994

Raw and derived transfer-function values¶

Use the numerical accessors for custom analyses. Magnitude describes gain, phase describes frequency-dependent alignment, and group delay summarizes the slope of unwrapped phase.

In [4]:

frequencies, transfer = model.transfer_function_at()
components = model.transfer_function_components_at(phase_unit="deg")

print("complex transfer values:", transfer.shape)
print("magnitude:", components.magnitude.shape)
print("phase unit:", components.phase_unit)
print("group-delay unit: seconds")

frequencies, transfer = model.transfer_function_at() components = model.transfer_function_components_at(phase_unit="deg") print("complex transfer values:", transfer.shape) print("magnitude:", components.magnitude.shape) print("phase unit:", components.phase_unit) print("group-delay unit: seconds")

complex transfer values: (65,)
magnitude: (65,)
phase unit: deg
group-delay unit: seconds

Transfer magnitude¶

In [5]:

model.plot_transfer_function(
    kind="magnitude",
    title="Gain of the learned mapping",
)
plt.show()

model.plot_transfer_function( kind="magnitude", title="Gain of the learned mapping", ) plt.show()

Transfer phase¶

In [6]:

model.plot_transfer_function(
    kind="phase",
    phase_unit="deg",
    title="Unwrapped transfer phase",
)
plt.show()

model.plot_transfer_function( kind="phase", phase_unit="deg", title="Unwrapped transfer phase", ) plt.show()

Group delay¶

Group delay is most interpretable over frequencies with appreciable transfer magnitude. Large values where gain is near zero are often numerically unstable and should not be overinterpreted.

In [7]:

model.plot_transfer_function(
    kind="group_delay",
    group_delay_unit="ms",
    title="Frequency-dependent group delay",
)
plt.show()

model.plot_transfer_function( kind="group_delay", group_delay_unit="ms", title="Frequency-dependent group delay", ) plt.show()

Compute held-out cross-spectral diagnostics¶

This object is reusable across coherence and cross-spectrum plots, avoiding repeated prediction and spectral estimation.

In [8]:

diagnostics = model.cross_spectral_diagnostics(
    stimulus=test_x,
    response=test_y,
)
print("coherence:", diagnostics.coherence.shape)
print("cross spectrum:", diagnostics.cross_spectrum.shape)

diagnostics = model.cross_spectral_diagnostics( stimulus=test_x, response=test_y, ) print("coherence:", diagnostics.coherence.shape) print("cross spectrum:", diagnostics.cross_spectrum.shape)

coherence: (65, 1)
cross spectrum: (65, 1)

Prediction-observation coherence¶

Coherence measures frequency-specific linear agreement. It is bounded between zero and one but does not test calibration.

In [9]:

model.plot_coherence(
    diagnostics=diagnostics,
    title="Held-out prediction coherence",
)
plt.show()

model.plot_coherence( diagnostics=diagnostics, title="Held-out prediction coherence", ) plt.show()

Predicted-observed cross spectrum¶

Magnitude shows frequency-specific covariation; phase shows their relative alignment. Plot them separately while interpreting the result because their units and scales differ.

In [10]:

model.plot_cross_spectrum(
    diagnostics=diagnostics,
    kind="magnitude",
    title="Held-out cross-spectrum magnitude",
)
plt.show()

model.plot_cross_spectrum( diagnostics=diagnostics, kind="magnitude", title="Held-out cross-spectrum magnitude", ) plt.show()

In [11]:

model.plot_cross_spectrum(
    diagnostics=diagnostics,
    kind="phase",
    phase_unit="deg",
    title="Held-out cross-spectrum phase",
)
plt.show()

model.plot_cross_spectrum( diagnostics=diagnostics, kind="phase", phase_unit="deg", title="Held-out cross-spectrum phase", ) plt.show()

Continue with the Diagnostics and Transfer Functions guide for a compact tool-selection summary.