Trial Weighting¶

Trial weighting changes how much each trial contributes to the fitted spectral statistics. It is appropriate only when the weight reflects a defensible, prespecified quality measure. High neural variance is not automatically noise.

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

from fftrf import inverse_variance_weights

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 from fftrf import inverse_variance_weights

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

Create a controlled case with one known noisy trial¶

Independent noise is injected deliberately, making response variance a valid noise proxy in this simulation.

In [3]:

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials(
    seed=77,
    n_trials=7,
    noise_scale=0.15,
)
train_x, test_x = stimulus[:-1], stimulus[-1]
train_y = [trial.copy() for trial in response[:-1]]
test_y = response[-1]

noise_rng = np.random.default_rng(0)
train_y[0] += 20.0 * noise_rng.standard_normal(train_y[0].shape)
weights = inverse_variance_weights(train_y)
print("trial weights:", np.round(weights, 3))

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials( seed=77, n_trials=7, noise_scale=0.15, ) train_x, test_x = stimulus[:-1], stimulus[-1] train_y = [trial.copy() for trial in response[:-1]] test_y = response[-1] noise_rng = np.random.default_rng(0) train_y[0] += 20.0 * noise_rng.standard_normal(train_y[0].shape) weights = inverse_variance_weights(train_y) print("trial weights:", np.round(weights, 3))

trial weights: [0.002 0.183 0.202 0.193 0.216 0.204]

Compare unweighted and weighted fits on the same held-out trial¶

In [4]:

common = dict(
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=1e-2,
    segment_length=512,
    overlap=0.5,
    window="hann",
)

unweighted = TRF(direction=1)
unweighted.train(train_x, train_y, **common)
_, unweighted_r = unweighted.predict(test_x, test_y)

weighted = TRF(direction=1)
weighted.train(train_x, train_y, trial_weights=weights, **common)
_, weighted_r = weighted.predict(test_x, test_y)

print(f"unweighted held-out r: {float(unweighted_r):.3f}")
print(f"weighted held-out r: {float(weighted_r):.3f}")

common = dict( fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_length=512, overlap=0.5, window="hann", ) unweighted = TRF(direction=1) unweighted.train(train_x, train_y, **common) _, unweighted_r = unweighted.predict(test_x, test_y) weighted = TRF(direction=1) weighted.train(train_x, train_y, trial_weights=weights, **common) _, weighted_r = weighted.predict(test_x, test_y) print(f"unweighted held-out r: {float(unweighted_r):.3f}") print(f"weighted held-out r: {float(weighted_r):.3f}")

unweighted held-out r: 0.974
weighted held-out r: 0.997

Inspect what the weighting changed¶

The controlled example should downweight trial 1 and move the recovered kernel toward the known signal. In real data, derive weights from an external artifact or noise measure rather than from whichever choice improves the test score.

In [5]:

fig, axes = plt.subplots(2, 1, figsize=(8, 6), constrained_layout=True)
axes[0].bar(np.arange(1, len(weights) + 1), weights)
axes[0].set(xlabel="Training trial", ylabel="Weight",
            title="Prespecified trial contributions")
axes[1].plot(times * 1e3, true_kernel, "--", color="black", label="True")
axes[1].plot(unweighted.times * 1e3, unweighted.weights[0, :, 0],
             label="Unweighted")
axes[1].plot(weighted.times * 1e3, weighted.weights[0, :, 0],
             label="Weighted")
axes[1].set(xlabel="Lag (ms)", ylabel="Weight", title="Recovered kernels")
axes[1].legend()
plt.show()

fig, axes = plt.subplots(2, 1, figsize=(8, 6), constrained_layout=True) axes[0].bar(np.arange(1, len(weights) + 1), weights) axes[0].set(xlabel="Training trial", ylabel="Weight", title="Prespecified trial contributions") axes[1].plot(times * 1e3, true_kernel, "--", color="black", label="True") axes[1].plot(unweighted.times * 1e3, unweighted.weights[0, :, 0], label="Unweighted") axes[1].plot(weighted.times * 1e3, weighted.weights[0, :, 0], label="Weighted") axes[1].set(xlabel="Lag (ms)", ylabel="Weight", title="Recovered kernels") axes[1].legend() plt.show()

See Trial Weighting and Bootstrap before applying this option to empirical neural data.