Decoding of stimuli from neural data

Here we will simulate neural data given a ground truth encoding model

and try to decode the stimulus from the data.

# Set up a neural model
from braincoder.models import GaussianPRF
import numpy as np
import pandas as pd
import scipy.stats as ss

# Set up 100 random of PRF parameters
n = 20
n_trials = 50
noise = 1.

mu = np.random.rand(n) * 100
sd = np.random.rand(n) * 45 + 5
amplitude = np.random.rand(n) * 5
baseline = np.random.rand(n) * 2 - 1

parameters = pd.DataFrame({'mu':mu, 'sd':sd, 'amplitude':amplitude, 'baseline':baseline})

# We have a paradigm of random numbers between 0 and 100
paradigm = np.ceil(np.random.rand(n_trials) * 100)

model = GaussianPRF(parameters=parameters)
data = model.simulate(paradigm=paradigm, noise=noise)

# Now we fit back the PRF parameters
from braincoder.optimize import ParameterFitter, ResidualFitter
fitter = ParameterFitter(model, data, paradigm)
mu_grid = np.arange(0, 100, 5)
sd_grid = np.arange(5, 50, 5)

grid_pars = fitter.fit_grid(mu_grid, sd_grid, [1.0], [0.0], use_correlation_cost=True, progressbar=False)
grid_pars = fitter.refine_baseline_and_amplitude(grid_pars)

for par in ['mu', 'sd', 'amplitude', 'baseline']:
    print(f'Correlation grid-fitted parameter and ground truth for *{par}*: {ss.pearsonr(grid_pars[par], parameters[par])[0]:0.2f}')

gd_pars = fitter.fit(init_pars=grid_pars, progressbar=False)

for par in ['mu', 'sd', 'amplitude', 'baseline']:
    print(f'Correlation gradient descent-fitted parameter and ground truth for *{par}*: {ss.pearsonr(grid_pars[par], parameters[par])[0]:0.2f}')

  0%|          | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:00<00:00, 302.58it/s]
Correlation grid-fitted parameter and ground truth for *mu*: 0.84
Correlation grid-fitted parameter and ground truth for *sd*: 0.09
Correlation grid-fitted parameter and ground truth for *amplitude*: 0.84
Correlation grid-fitted parameter and ground truth for *baseline*: 0.52
Correlation gradient descent-fitted parameter and ground truth for *mu*: 0.84
Correlation gradient descent-fitted parameter and ground truth for *sd*: 0.09
Correlation gradient descent-fitted parameter and ground truth for *amplitude*: 0.84
Correlation gradient descent-fitted parameter and ground truth for *baseline*: 0.52

# Now we fit the covariance matrix
stimulus_range = np.arange(1, 100).astype(np.float32)

model.init_pseudoWWT(stimulus_range=stimulus_range, parameters=gd_pars)
resid_fitter = ResidualFitter(model, data, paradigm, gd_pars)
omega, dof = resid_fitter.fit(progressbar=False)

init_tau: 0.7906791567802429, 1.191803216934204
USING A PSEUDO-WWT!
WWT max: 1339.4237060546875

# Now we simulate unseen test data:
test_paradigm = np.ceil(np.random.rand(n_trials) * 100)
test_data = model.simulate(paradigm=test_paradigm, noise=noise)

# And decode the test paradigm
posterior = model.get_stimulus_pdf(test_data, stimulus_range, model.parameters, omega=omega, dof=dof)
# Rename the row index so it doesn't clash with the column index (both default to 'stimulus')
posterior.index.name = 'frame'

# Finally, we make some plots to see how well the decoder did
import matplotlib.pyplot as plt
import seaborn as sns

tmp = posterior.set_index(pd.Series(test_paradigm, name='ground truth'), append=True).loc[:8].stack().to_frame('p')

g = sns.FacetGrid(tmp.reset_index(), col='frame', col_wrap=3)

g.map(plt.plot, 'stimulus', 'p', color='k')

def test(data, **kwargs):
    plt.axvline(data.mean(), c='k', ls='--', **kwargs)
g.map(test, 'ground truth')
g.set(xlabel='Stimulus value', ylabel='Posterior probability density')

<seaborn.axisgrid.FacetGrid object at 0x30ab36510>

# Let's look at the summary statistics of the posteriors posteriors
def get_posterior_stats(posterior, normalize=True):
    posterior = posterior.copy()
    posterior = posterior.div(np.trapezoid(posterior, posterior.columns,axis=1), axis=0)

    # Take integral over the posterior to get to the expectation (mean posterior)
    E = np.trapezoid(posterior*posterior.columns.values[np.newaxis,:], posterior.columns, axis=1)

    # Take the integral over the posterior to get the expectation of the distance to the
    # mean posterior (i.e., standard deviation)
    sd = np.trapezoid(np.abs(E[:, np.newaxis] - posterior.columns.astype(float).values[np.newaxis, :]) * posterior, posterior.columns, axis=1)

    stats = pd.DataFrame({'E':E, 'sd':sd}, index=posterior.index)
    return stats

posterior_stats = get_posterior_stats(posterior)

# Let's see how far the posterior mean is from the ground truth
plt.errorbar(test_paradigm, posterior_stats['E'],posterior_stats['sd'], fmt='o',)
plt.plot([0, 100], [0,100], c='k', ls='--')

plt.xlabel('Ground truth')
plt.ylabel('Mean posterior')

# Let's see how the error depends on the standard deviation of the posterior
error = test_paradigm - posterior_stats['E']
error_abs = np.abs(error)
error_abs.name = 'error'

sns.lmplot(x='sd', y='error', data=posterior_stats.join(error_abs))

plt.xlabel('Standard deviation of posterior')
plt.ylabel('Objective error')

• 
• 

Text(29.0, 0.5, 'Objective error')

# Now, let's try to find the MAP estimate using gradient descent
from braincoder.optimize import StimulusFitter
stimulus_fitter = StimulusFitter(model=model, data=test_data, omega=omega)

# We start with a very coarse grid search, so we are sure we are in the right ballpark
estimated_stimuli_grid = stimulus_fitter.fit_grid(np.arange(1, 100, 5))

# We can then refine the estimate using gradient descent
estimated_stimuli_gd = stimulus_fitter.fit(init_pars=estimated_stimuli_grid, progressbar=False)

# Let's see how well we did
plt.scatter(test_paradigm, estimated_stimuli_grid, alpha=.5, label='MAP (grid search)')
plt.scatter(test_paradigm, estimated_stimuli_gd, alpha=.5, label='MAP (gradient descent)')
plt.scatter(test_paradigm, posterior_stats['E'], alpha=.5, label='Mean posterior')
plt.plot([0, 100], [0,100], c='k', ls='--', label='Identity line')
plt.legend()
# %%

  0%|          | 0/1000 [00:00<?, ?it/s]
  0%|          | 1/1000 [00:00<01:45,  9.44it/s]
  3%|▎         | 31/1000 [00:00<00:05, 176.27it/s]
  6%|▋         | 63/1000 [00:00<00:03, 239.63it/s]
 10%|▉         | 97/1000 [00:00<00:03, 276.20it/s]
 13%|█▎        | 130/1000 [00:00<00:02, 293.87it/s]
 16%|█▌        | 162/1000 [00:00<00:02, 301.37it/s]
 20%|█▉        | 195/1000 [00:00<00:02, 308.01it/s]
 23%|██▎       | 228/1000 [00:00<00:02, 313.27it/s]
 26%|██▌       | 262/1000 [00:00<00:02, 320.96it/s]
 30%|██▉       | 297/1000 [00:01<00:02, 329.37it/s]
 33%|███▎      | 333/1000 [00:01<00:01, 337.68it/s]
 37%|███▋      | 369/1000 [00:01<00:01, 342.09it/s]
 40%|████      | 405/1000 [00:01<00:01, 346.12it/s]
 44%|████▍     | 441/1000 [00:01<00:01, 349.33it/s]
 48%|████▊     | 476/1000 [00:01<00:01, 347.72it/s]
 51%|█████     | 511/1000 [00:01<00:01, 345.18it/s]
 55%|█████▍    | 546/1000 [00:01<00:01, 338.45it/s]
 58%|█████▊    | 581/1000 [00:01<00:01, 339.05it/s]
 62%|██████▏   | 617/1000 [00:01<00:01, 343.93it/s]
 65%|██████▌   | 652/1000 [00:02<00:01, 337.82it/s]
 69%|██████▊   | 687/1000 [00:02<00:00, 339.83it/s]
 72%|███████▏  | 723/1000 [00:02<00:00, 343.48it/s]
 76%|███████▌  | 758/1000 [00:02<00:00, 342.07it/s]
 79%|███████▉  | 793/1000 [00:02<00:00, 344.25it/s]
 83%|████████▎ | 828/1000 [00:02<00:00, 344.81it/s]
 86%|████████▋ | 863/1000 [00:02<00:00, 345.62it/s]
 90%|████████▉ | 898/1000 [00:02<00:00, 339.64it/s]
 93%|█████████▎| 933/1000 [00:02<00:00, 340.71it/s]
 97%|█████████▋| 968/1000 [00:02<00:00, 340.78it/s]
100%|██████████| 1000/1000 [00:03<00:00, 326.34it/s]

<matplotlib.legend.Legend object at 0x30e86c110>