Getting Started Notebook¶

This notebook walks through a small but meaningful forward-model example: a simulated speech-envelope regressor predicts one simulated brain-response channel. The data are synthetic, but the held-out score and recovered kernel are informative rather than random.

In [1]:

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

from fftrf import TRF, inverse_variance_weights


def lag_times(fs: float, tmin: float, tmax: float) -> np.ndarray:
    lag_start = int(round(tmin * fs))
    lag_stop = int(round(tmax * fs))
    return np.arange(lag_start, lag_stop, dtype=int) / fs


def make_speech_envelope(
    n_samples: int,
    fs: float,
    rng: np.random.Generator,
) -> np.ndarray:
    cutoff = min(8.0 / (0.5 * fs), 0.99)
    b, a = butter(3, cutoff, btype="low")
    envelope = filtfilt(b, a, rng.standard_normal(n_samples))
    envelope = envelope - envelope.min()
    envelope = envelope / np.clip(envelope.std(), np.finfo(float).eps, None)
    return envelope


def simulate_trial(
    rng: np.random.Generator,
    *,
    n_samples: int,
    fs: float,
    kernel: np.ndarray,
    noise_scale: float,
) -> tuple[np.ndarray, np.ndarray]:
    envelope = make_speech_envelope(n_samples=n_samples, fs=fs, rng=rng)
    response = np.convolve(envelope, kernel, mode="full")[:n_samples]
    response += noise_scale * rng.standard_normal(n_samples)
    return envelope[:, np.newaxis], response[:, np.newaxis]


rng = np.random.default_rng(7)
fs = 128.0
tmin = 0.0
tmax = 0.320
times = lag_times(fs, tmin, tmax)
true_kernel = (
    0.90 * np.exp(-0.5 * ((times - 0.050) / 0.015) ** 2)
    - 0.55 * np.exp(-0.5 * ((times - 0.115) / 0.022) ** 2)
    + 0.18 * np.exp(-0.5 * ((times - 0.210) / 0.035) ** 2)
)

trials = [
    simulate_trial(
        rng,
        n_samples=3_072,
        fs=fs,
        kernel=true_kernel,
        noise_scale=0.10 + 0.03 * rng.random(),
    )
    for _ in range(6)
]
stimulus = [trial_stimulus for trial_stimulus, _ in trials]
response = [trial_response for _, trial_response in trials]

train_stimulus = stimulus[:-1]
train_response = response[:-1]
test_stimulus = stimulus[-1]
test_response = response[-1]

model = TRF(direction=1)
cv_scores = model.train(
    stimulus=train_stimulus,
    response=train_response,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=np.logspace(-6, 0, 7),
    segment_duration=2.0,
    overlap=0.5,
    window="hann",
    k="loo",
    trial_weights=inverse_variance_weights(train_response),
)

_, score = model.predict(stimulus=test_stimulus, response=test_response)
kernel_correlation = np.corrcoef(true_kernel, model.weights[0, :, 0])[0, 1]

print(f"Selected regularization: {model.regularization}")
print(f"Held-out Pearson r: {float(score):.3f}")
print(f"Kernel correlation to ground truth: {kernel_correlation:.3f}")

import matplotlib matplotlib.use("module://matplotlib_inline.backend_inline") import matplotlib.pyplot as plt import numpy as np from scipy.signal import butter, filtfilt from fftrf import TRF, inverse_variance_weights def lag_times(fs: float, tmin: float, tmax: float) -> np.ndarray: lag_start = int(round(tmin * fs)) lag_stop = int(round(tmax * fs)) return np.arange(lag_start, lag_stop, dtype=int) / fs def make_speech_envelope( n_samples: int, fs: float, rng: np.random.Generator, ) -> np.ndarray: cutoff = min(8.0 / (0.5 * fs), 0.99) b, a = butter(3, cutoff, btype="low") envelope = filtfilt(b, a, rng.standard_normal(n_samples)) envelope = envelope - envelope.min() envelope = envelope / np.clip(envelope.std(), np.finfo(float).eps, None) return envelope def simulate_trial( rng: np.random.Generator, *, n_samples: int, fs: float, kernel: np.ndarray, noise_scale: float, ) -> tuple[np.ndarray, np.ndarray]: envelope = make_speech_envelope(n_samples=n_samples, fs=fs, rng=rng) response = np.convolve(envelope, kernel, mode="full")[:n_samples] response += noise_scale * rng.standard_normal(n_samples) return envelope[:, np.newaxis], response[:, np.newaxis] rng = np.random.default_rng(7) fs = 128.0 tmin = 0.0 tmax = 0.320 times = lag_times(fs, tmin, tmax) true_kernel = ( 0.90 * np.exp(-0.5 * ((times - 0.050) / 0.015) ** 2) - 0.55 * np.exp(-0.5 * ((times - 0.115) / 0.022) ** 2) + 0.18 * np.exp(-0.5 * ((times - 0.210) / 0.035) ** 2) ) trials = [ simulate_trial( rng, n_samples=3_072, fs=fs, kernel=true_kernel, noise_scale=0.10 + 0.03 * rng.random(), ) for _ in range(6) ] stimulus = [trial_stimulus for trial_stimulus, _ in trials] response = [trial_response for _, trial_response in trials] train_stimulus = stimulus[:-1] train_response = response[:-1] test_stimulus = stimulus[-1] test_response = response[-1] model = TRF(direction=1) cv_scores = model.train( stimulus=train_stimulus, response=train_response, fs=fs, tmin=tmin, tmax=tmax, regularization=np.logspace(-6, 0, 7), segment_duration=2.0, overlap=0.5, window="hann", k="loo", trial_weights=inverse_variance_weights(train_response), ) _, score = model.predict(stimulus=test_stimulus, response=test_response) kernel_correlation = np.corrcoef(true_kernel, model.weights[0, :, 0])[0, 1] print(f"Selected regularization: {model.regularization}") print(f"Held-out Pearson r: {float(score):.3f}") print(f"Kernel correlation to ground truth: {kernel_correlation:.3f}")

Selected regularization: 0.01
Held-out Pearson r: 0.999
Kernel correlation to ground truth: 0.997

What to inspect first¶

• model.weights is the recovered lag-domain TRF.
• model.regularization is the lambda selected by leave-one-out cross-validation.
• score is the held-out correlation on a trial the model never saw during fitting.

In [2]:

diagnostics = model.cross_spectral_diagnostics(
    stimulus=test_stimulus,
    response=test_response,
    trial_weights=None,
)

fig, axes = plt.subplots(2, 1, figsize=(9, 7), constrained_layout=True)
model.plot(
    input_index=0,
    output_index=0,
    ax=axes[0],
    title=f"Recovered speech-envelope TRF (held-out r = {float(score):.3f})",
)
model.plot_coherence(
    diagnostics=diagnostics,
    output_index=0,
    ax=axes[1],
    title="Held-out prediction coherence",
)
plt.show()

diagnostics = model.cross_spectral_diagnostics( stimulus=test_stimulus, response=test_response, trial_weights=None, ) fig, axes = plt.subplots(2, 1, figsize=(9, 7), constrained_layout=True) model.plot( input_index=0, output_index=0, ax=axes[0], title=f"Recovered speech-envelope TRF (held-out r = {float(score):.3f})", ) model.plot_coherence( diagnostics=diagnostics, output_index=0, ax=axes[1], title="Held-out prediction coherence", ) plt.show()

Related next steps¶

• Read the full Getting Started page for parameter-by-parameter guidance.
• Walk through the intended lifecycle in Core Workflow.
• Open Examples to find the script that most closely matches your data.
• Use score(...), plot_grid(...), cross_spectral_diagnostics(...), or bootstrap_confidence_interval(...) once the first fit is working.