Optimizing hyperparameters using L-BFGS-B¶

# Required imports
import numpy as np
from scipy.stats import zscore, norm
import mne
from posthoc import (Workbench, WorkbenchOptimizer, cov_estimators,
                     normalizers)
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import accuracy_score, roc_auc_score
from sklearn.linear_model import LogisticRegression
from matplotlib import pyplot as plt

import warnings
warnings.simplefilter('ignore')

path = mne.datasets.sample.data_path()
raw = mne.io.read_raw_fif(path + '/MEG/sample/sample_audvis_raw.fif',
                          preload=True)
events = mne.find_events(raw)
event_id = dict(left=1, right=2)
raw.pick_types(meg='grad')
raw.filter(None, 20)
raw, events = raw.resample(50, events=events)

# Create epochs
epochs = mne.Epochs(raw, events, event_id, tmin=-0.2, tmax=0.5,
                    baseline=(-0.2, 0), preload=True)
n_epochs, n_channels, n_samples = epochs.get_data().shape

X = epochs.get_data().reshape(len(epochs), -1)

# There is currently a bug in the intercept calculation for the
# kernel-formulated leave-one-out code path (for the paper, the slower
# 'traditional' path was used). For this example, we'll just zscore the data
# and abandon the intercept calculation altogether.
X = zscore(X, axis=0)

# Create training labels, based on the event codes during the experiment.
y = epochs.events[:, [2]]
y = y - 1.5
y *= 2
y = y.astype(int)

# Split the data 50/50, but make sure the number of left/right epochs are
# balanced.
folds = StratifiedKFold(n_splits=2)
train_index, test_index = next(folds.split(X, y))
X_train, y_train = X[train_index], y[train_index]
X_test, y_test = X[test_index], y[test_index]

# The logistic regression model ignores observations that are close to the
# decision boundary. The parameter `C` controls how far away observations have
# to be in order to not be ignored. A setting of 25 means "quite far". We also
# specify the seed for the random number generator, so that this example
# replicates exactly every time.
base_model = LogisticRegression(C=25, solver='lbfgs', random_state=0,
                                fit_intercept=False)

# Train on the training data and predict the test data.
base_model.fit(X_train, y_train)
y_hat = base_model.predict(X_test)

# How many epochs did we decode correctly?
base_model_accuracy = accuracy_score(y_test, y_hat)
print('Base model accuracy:', base_model_accuracy)

base_model = Workbench(base_model).fit(X_train, y_train)

# Plot the pattern
plt.figure()
plt.plot(epochs.times, base_model.pattern_.reshape(n_channels, n_samples).T,
         color='black', alpha=0.2)
plt.xlabel('Time (s)')
plt.ylabel('Signal (normalized units)')
plt.title('Pattern learned by the base model')

# The function that modifies the pattern takes as input the original pattern,
# the training data, and two parameters that define the center and width of the
# Gaussian kernel.
cache = dict()
def pattern_modifier(pattern, X, y, center, width):
    """Multiply the pattern with a Gaussian kernel."""
    mod_pattern = pattern.reshape(n_channels, n_samples)
    key = (center, width)
    if key in cache:
        kernel = cache[key]
    else:
        kernel = norm(center, width).pdf(np.arange(n_samples))
        kernel /= kernel.max()
        cache[key] = kernel
    mod_pattern = mod_pattern * kernel[np.newaxis, :]
    return mod_pattern.reshape(pattern.shape)

def scorer(model, X, y):
    """Our scoring function."""
    y_hat = model.predict(X)
    return roc_auc_score(y, y_hat)

# Define initial values for the parameters. The optimization algorithm will
# use gradient descend using these as starting point.
initial_center = np.searchsorted(epochs.times, 0.05)
initial_width = 5

# Define the allowed range for the parameters. The optimizer will not exceed
# these.
center_bounds = (5, 25)
width_bounds = (1, 50)

# Define the post-hoc model using an optimizer to fine-tune the parameters.
optimized_model = WorkbenchOptimizer(
    base_model,
    cov=cov_estimators.ShrinkageKernel(1.0),
    cov_param_bounds=[(0.9, 1.0)],
    pattern_modifier=pattern_modifier,
    pattern_param_x0=[initial_center, initial_width],
    pattern_param_bounds=[center_bounds, width_bounds],
    normalizer_modifier=normalizers.unit_gain,
    scoring=scorer,
    verbose=True,
    random_search=20,
).fit(X_train, y_train)

# Decode the test data
y_hat = optimized_model.predict(X_test).ravel()

# Assign the 'left' class to values above 0 and 'right' to values below 0
y_bin = np.zeros(len(y_hat), dtype=np.int)
y_bin[y_hat >= 0] = 1
y_bin[y_hat < 0] = -1

# How many epochs did we decode correctly?
optimized_model_accuracy = accuracy_score(y_test, y_bin)
print('Base model accuracy:', base_model_accuracy)
print('Optimized model accuracy:', optimized_model_accuracy)

plt.figure()
plt.plot(epochs.times,
         optimized_model.pattern_.reshape(n_channels, n_samples).T,
         color='black', alpha=0.2)
plt.xlabel('Time (s)')
plt.ylabel('Signal (normalized units)')
plt.title('Optimized pattern')