Reproducing Makin et al. (2018) from dandiset 000129: decoders at 64 ms¶
Paper. Makin JG, O'Doherty JE, Cardoso MMB, Sabes PN (2018). Superior arm-movement decoding from cortex with a new, unsupervised-learning algorithm. J. Neural Eng. 15, 026010. doi:10.1088/1741-2552/aa9e95
Dataset. DANDI:000129 version 0.241017.1444 (doi:10.48324/dandi.000129/0.241017.1444), "MC_RTT", streamed from the archive. The paper cites no version of the dandiset; 0.241017.1444 is the only published version (the draft holds the same two assets, 50,965,512 bytes, and differs only in its title), so it is pinned.
Claims. claims/000129/10.1088-1741-2552-aa9e95.yaml, written from the version of record
(main text, Methods and legends; the paper has no supplement) and the paper's code
(github.com/jgmakin/rbmish at commit 700062f5).
What the dataset holds. One session of monkey I (Indy), M1 array only, from the Neural Latents Benchmark: a training file with 649.1 s of sorted-unit spikes and kinematics, and a test file with 190.3 s of spikes whose kinematics are withheld. The paper used 49 sessions from three monkeys. Claims judged on a mean over sessions are tested here on this one session (a deviation); claims that need the bootstrap across the 49 sessions, or the per-session spread, cannot be tested.
This notebook covers the decoders at 64 ms bins with 320 s of training (figures 2, 4 and 8).
The notebook binwidth_training_neurons.ipynb covers bin widths, amounts of training data and
fractions of neurons (figures 5 to 7).
| Claim | Here |
|---|---|
| units-monkey-i | tested |
| fig2-velocity-boost, fig2-kflatent-dynamic-little-effect, fig2-refh-static-vel-acc, fig2-position-dynamic-gain, fig2-sagittal-better, fig2-linear-worst, fig2-gross-failures | tested on one session |
| fig4-vel-acc-every-session (per-session part), fig4-best-variable, kflatent-modest | tested on one session |
| fig8-velocity-variants, fig8-gaussian-still-beats-ukf | tested on one session |
| fig5-advantage-with-binwidth, fig5-derivative-order, fig7-affine-position, fig7-kinematic-model-saturates, advantage-grows-with-derivative | tested in binwidth_training_neurons.ipynb |
| fig2-refh-kf-best, intro-gain-over-nearest, fig4-position, fig8-position-variants, static-refh-position-worse | blocked: bootstrap across the paper's 49 sessions, which the dandiset does not hold; the single-session values are printed below |
| fig4-advantage-grows, fig5-vel-acc-all-bins, fig5-position-bins, fig6-training-time, fig7-vel-acc-neurons, fig7-position-neurons, fig9-taps | blocked: bootstrap or regression across the paper's 49 sessions (single-session values of the fig5 to fig7 claims in binwidth_training_neurons.ipynb; the tapped decoders of fig9-taps are not implemented) |
| fig4-task-similar | blocked: monkey C's data and monkey I's other sessions |
import os
os.environ.setdefault("OMP_NUM_THREADS", "4")
import time
import numpy as np
import pandas as pd
import requests
import makin as M
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
T0 = time.time()
d = M.load_session()
print(f"dandiset {M.DANDISET} version {M.VERSION}")
print(d["session_description"])
print("session_start_time", d["session_start_time"])
dandiset 000129 version 0.241017.1444 Data from monkey Indy performing self-paced random target reaching task. This file contains continuous segments of the full session on 2017-02-02 that can be used for training models for the Neural Latents Benchmark. session_start_time 2017-02-02T00:00:00-08:00
Data check¶
Each check of the metadata against the values is stated under its output.
ids = d["unit_ids"]
n_spk = np.array([len(s) for s in d["spikes"]])
dur = len(d["finger_pos"]) / d["rate"]
print(f"subjects 1 (Indy), sessions 1, units {len(ids)} ({d['heldout'].sum()} marked heldout by the benchmark)")
print(f"kinematics: {len(d['finger_pos'])} samples at {d['rate']:.0f} Hz from t0 = {d['t0']} s, duration {dur:.1f} s")
print("finger_pos unit attribute:", d["unit"], " conversion:", d["conversion"])
for k in ["finger_pos", "cursor_pos", "target_pos"]:
a = d[k]
print(f"{k:11s} x {np.nanmin(a[:,0]):7.1f} to {np.nanmax(a[:,0]):6.1f}, y {np.nanmin(a[:,1]):6.1f} to {np.nanmax(a[:,1]):6.1f} (stored values); NaN samples {np.isnan(a[:,0]).sum()}")
gaps = M.segments(np.isnan(d["finger_pos"][:, 0]))
print("NaN gaps (s):", [(round(s / d['rate'], 3), round(e / d['rate'], 3)) for s, e in gaps])
st_all = np.concatenate(d["spikes"])
print(f"spikes {len(st_all)}, times {st_all.min():.3f} to {st_all.max():.3f} s, sorted within units: {all(np.all(np.diff(s) >= 0) for s in d['spikes'])}")
in_gap = sum(((st_all >= s / d['rate']) & (st_all < e / d['rate'])).sum() for s, e in gaps)
print("spikes inside the NaN gaps:", in_gap)
print("electrode locations:", pd.Series(d["electrode_locations"]).value_counts().to_dict())
print("successive 1 kHz finger samples that are equal:", int((np.diff(d["finger_pos"][:, 0]) == 0).sum()))
print("target x positions (stored):", np.unique(d["target_pos"][~np.isnan(d["target_pos"][:, 0]), 0]))
print("unit ids end in:", sorted(set(int(i) % 100 for i in ids)), "(channel*100 + unit; 0 would be the unsorted hash unit)")
print("channels with units:", len(set(int(i) // 100 for i in ids)))
rate = n_spk / dur
print(f"rates (Hz) min {rate.min():.2f}, median {np.median(rate):.2f}, max {rate.max():.2f}")
print("trials table splits:", pd.Series(d["trials_split"]).value_counts().to_dict())
subjects 1 (Indy), sessions 1, units 130 (32 marked heldout by the benchmark)
kinematics: 649100 samples at 1000 Hz from t0 = 0.0 s, duration 649.1 s
finger_pos unit attribute: meters conversion: 0.001
finger_pos x -63.8 to 70.6, y -7.5 to 117.8 (stored values); NaN samples 300
cursor_pos x -65.2 to 71.0, y -7.8 to 117.9 (stored values); NaN samples 300
target_pos x -52.5 to 52.5, y 7.5 to 112.5 (stored values); NaN samples 300
NaN gaps (s): [(np.float64(162.2), np.float64(162.3)), (np.float64(324.5), np.float64(324.6)), (np.float64(486.8), np.float64(486.9))]
spikes 339874, times 0.002 to 649.096 s, sorted within units: True
spikes inside the NaN gaps: 0
electrode locations: {'M1': 96}
successive 1 kHz finger samples that are equal: 0
target x positions (stored): [-52.5 -37.5 -22.5 -7.5 7.5 22.5 37.5 52.5]
unit ids end in: [1, 2, 3, 4] (channel*100 + unit; 0 would be the unsorted hash unit)
channels with units: 82
rates (Hz) min 0.08, median 2.96, max 19.51
trials table splits: {'train': 810, 'val': 270}
The kinematic series are stored in mm with unit "meters" and conversion 0.001, which is consistent (value times conversion gives meters); the values span about 135 by 125 mm, and the targets sit at eight x positions 15 mm apart, as the 8 by 8 grid of the description and of the paper's later monkey I sessions has. All electrodes are labelled M1. The finger position description says 250 Hz upsampled to 1 kHz; successive 1 kHz samples are never equal, so the upsampling interpolated, and taking every fourth sample recovers a 250 Hz series. Three gaps of 0.1 s are NaN in every behavioural series. The unit ids end in 1 to 4: no unsorted hash unit is in the file, though the paper kept one per channel (Methods, Data acquisition).
# Is this one of the paper's sessions? The paper's session list (code, NeuralAnalysis/Indy_HC_Sessions.csv)
CSV = ("https://raw.githubusercontent.com/jgmakin/rbmish/700062f5118b0910411cf9a8f73498e22203731d/"
"NeuralAnalysis/Indy_HC_Sessions.csv")
sess = pd.read_csv(CSV, header=None, names=["date", "block", "n_reaches", "arrays", "grid", "use"])
print(f"{len(sess)} rows; dates {sess.date.min()} to {sess.date.max()}; rows with at least 300 reaches: {(sess.n_reaches >= 300).sum()}")
print("2017-02-02 in the list:", 20170202 in set(sess.date))
print(sess.tail(4).to_string(index=False))
48 rows; dates 20160407 to 20170131; rows with at least 300 reaches: 37
2017-02-02 in the list: False
date block n_reaches arrays grid use
20170123 2 500 M1 8x8 y
20170124 1 500 M1 8x8 y
20170127 3 600 M1 8x8 y
20170131 2 650 M1 8x8 y
The file dates the session 2017-02-02. The paper's session list ends on 2017-01-31 and does not hold that date, so the session in the dandiset is probably not one of the paper's 37 monkey I sessions (it could also be a misdated one; nothing in the file links it to the list). It matches the paper's description of monkey I's later sessions: M1 array only, 8 by 8 grid.
# Test file: spikes only (kinematics withheld), and whether it repeats any of the training file
ft = M.open_asset(M.TEST_ASSET)
tid = ft["units/id"][:]
tst = np.split(ft["units/spike_times"][:], ft["units/spike_times_index"][:-1])
print("test file units:", len(tid), " behaviour in test file:", "processing" in ft and len(ft["processing"]) > 0)
print(f"test file spikes {sum(len(s) for s in tst)}, span {min(s.min() for s in tst if len(s)):.2f} to {max(s.max() for s in tst if len(s)):.2f} s")
# duplicated recording check: runs of 4 inter-spike intervals (to 1 us) shared between the files, per unit
idmap = {int(i): k for k, i in enumerate(ids)}
shared = total = 0
for i, s in zip(tid, tst):
a = d["spikes"][idmap[int(i)]]
keys = lambda x: set(map(tuple, np.round(np.lib.stride_tricks.sliding_window_view(np.diff(x), 4) * 1e6).astype(np.int64)))
kt = keys(s); total += len(kt); shared += len(kt & keys(a))
print(f"4-ISI patterns in the test file {total}, also in the training file {shared}")
ft.close()
test file units: 98 behaviour in test file: False test file spikes 66548, span 0.00 to 190.30 s
4-ISI patterns in the test file 66156, also in the training file 0
The test file holds 98 of the 130 units (the benchmark's held-in units) and no kinematics, so it cannot be scored; only the training file is used. No run of four inter-spike intervals is shared between the files, so the test file does not repeat the training file.
# Clock check: decode velocity from spike counts with the kinematics shifted in time.
R, X, edges = M.bin_data(d, 64)
keep = M.unit_filter(R, 64, 5000)
ok = ~np.isnan(X).any(1)
rows = []
for lag in range(-8, 9): # kinematics at bin t+lag decoded from spikes at bin t
Xs = np.roll(X, -lag, axis=0)
okl = ok & np.roll(ok, -lag) & (np.arange(len(X)) >= 8) & (np.arange(len(X)) < len(X) - 8)
tr = okl & (np.arange(len(X)) < 5000); te = okl & (np.arange(len(X)) >= 5000)
B = M.ridge(R[tr][:, keep], Xs[tr])
r2 = M.rsq(M.apply(B, R[te][:, keep]), Xs[te])
rows.append(dict(lag_ms=lag * 64, R2_vx=r2[2], R2_vy=r2[3], R2_x=r2[0]))
lagdf = pd.DataFrame(rows)
print(lagdf.round(3).to_string(index=False))
best = lagdf.loc[lagdf[["R2_vx", "R2_vy"]].mean(axis=1).idxmax(), "lag_ms"]
print("best lag for velocity (ms):", best)
lag_ms R2_vx R2_vy R2_x
-512 0.022 0.073 0.027
-448 0.023 0.078 0.032
-384 0.027 0.077 0.038
-320 0.038 0.084 0.044
-256 0.056 0.102 0.051
-192 0.079 0.127 0.059
-128 0.122 0.160 0.068
-64 0.195 0.206 0.079
0 0.274 0.275 0.096
64 0.328 0.327 0.122
128 0.312 0.320 0.157
192 0.217 0.268 0.194
256 0.108 0.182 0.223
320 0.045 0.099 0.237
384 0.020 0.054 0.237
448 0.017 0.040 0.227
512 0.019 0.036 0.213
best lag for velocity (ms): 64
Velocity is decoded best when the kinematics lag the spikes by about one to two bins, and the fit falls off on both sides; the spikes lead the movement, as motor cortex does, so the spike times and the kinematics share one clock. The paper's decoders use no offset (Methods, KFobs), and neither does this attempt.
print(f"bins of 64 ms: {len(X)}; training bins (first 320 s) 5000; test bins {len(X) - 5000}")
print(f"bins with NaN kinematics: {(~ok).sum()} (training {(~ok[:5000]).sum()}, test {(~ok[5000:]).sum()})")
rate64 = R.mean(0) / 0.064
dup = (rate64 > 0.5).sum() - keep.sum()
print(f"units: {R.shape[1]} in file, {(rate64 > 0.5).sum()} above 0.5 Hz, {dup} dropped as duplicates, {keep.sum()} used")
Xk = X[ok]
for i, v in enumerate(M.VARS):
print(f"{v:2s}: {Xk[:, i].min():9.1f} to {Xk[:, i].max():8.1f} (mm, mm/s, mm/s^2); test SD {np.nanstd(X[5000:, i]):.1f}")
bins of 64 ms: 10142; training bins (first 320 s) 5000; test bins 5142 bins with NaN kinematics: 6 (training 2, test 4) units: 130 in file, 111 above 0.5 Hz, 0 dropped as duplicates, 111 used x : -63.2 to 70.6 (mm, mm/s, mm/s^2); test SD 28.3 y : -7.4 to 117.7 (mm, mm/s, mm/s^2); test SD 27.6 vx: -466.1 to 463.9 (mm, mm/s, mm/s^2); test SD 73.2 vy: -324.7 to 346.8 (mm, mm/s, mm/s^2); test SD 67.4 ax: -3987.8 to 4977.3 (mm, mm/s, mm/s^2); test SD 685.0 ay: -3038.9 to 4348.0 (mm, mm/s, mm/s^2); test SD 547.5
Decoders at 64 ms¶
All eight decoders of the paper are fitted on the first 320 s and tested on the remaining 329 s
(the Methods, "Training and testing data"). The decoders and where each detail comes from are in
the docstrings of makin.py. The kinematics are the Methods reading (fingertip, 10 Hz
Butterworth, last sample of each bin); the code's reading (cursor, no filter, bin means) is run
as a variant below. Scores: R2 per variable on the test bins, and SNR = -10 log10(1 - R2) in dB.
def fit_config(seed=0, **kw):
sp = M.prepare(d, **kw)
t = time.time()
info = {}
pred = M.run_all(sp, seed, info=info)
r2 = pd.DataFrame({k: sp.score(v) for k, v in pred.items()}, index=M.VARS).T.loc[M.DECODERS]
return sp, pred, r2, time.time() - t, info
sp, pred, r2, secs, info = fit_config(seed=0)
print(f"fitted in {secs:.0f} s; units {sp.Rtr.shape[1]}, training bins {len(sp.Rtr)} ({sp.oktr.sum()} scored), test bins {len(sp.Rte)} ({sp.okte.sum()} scored)")
snr = M.snr(r2)
print("R2:"); print(r2.round(3).to_string())
print("\nSNR (dB):"); print(snr.round(2).to_string())
kfi = info["kflatent"]
print(f"\nKFlatent EM: {kfi['k']} latent states, {kfi['epochs']} E steps, log-likelihood {kfi['ll_first']:.0f} to {kfi['ll_last']:.0f}; "
f"gain of the last epoch {kfi['last_rel_gain']:.2e} of the total (criterion 1e-4); met the criterion: {bool(kfi['converged'])}")
print("static rEFH map, R2 on the training bins it was fitted on:", dict(zip(M.VARS, info["refh_train_r2"].round(3))))
print("static rEFH map, R2 on the test bins: ", dict(zip(M.VARS, r2.loc["refh"].values.round(3))))
fitted in 59 s; units 111, training bins 5000 (4998 scored), test bins 5142 (5138 scored)
R2:
x y vx vy ax ay
linear 0.096 0.234 0.274 0.275 0.055 0.116
kfobs 0.641 0.722 0.504 0.560 0.154 0.260
kflatent 0.561 0.651 0.502 0.538 0.150 0.251
kflatent_kf 0.579 0.666 0.528 0.557 0.166 0.269
wf 0.470 0.518 0.508 0.545 0.217 0.309
ukf 0.678 0.767 0.544 0.622 0.165 0.282
refh 0.189 0.380 0.594 0.598 0.326 0.370
refh_kf 0.703 0.778 0.680 0.709 0.397 0.444
SNR (dB):
x y vx vy ax ay
linear 0.44 1.16 1.39 1.40 0.25 0.53
kfobs 4.45 5.56 3.05 3.56 0.73 1.31
kflatent 3.58 4.58 3.03 3.36 0.71 1.26
kflatent_kf 3.75 4.76 3.26 3.54 0.79 1.36
wf 2.75 3.17 3.08 3.42 1.06 1.61
ukf 4.93 6.33 3.41 4.23 0.78 1.44
refh 0.91 2.08 3.92 3.96 1.71 2.01
refh_kf 5.27 6.53 4.95 5.37 2.19 2.55
KFlatent EM: 37 latent states, 100 E steps, log-likelihood -347650 to -338758; gain of the last epoch 4.25e-04 of the total (criterion 1e-4); met the criterion: False
static rEFH map, R2 on the training bins it was fitted on: {'x': np.float64(0.386), 'y': np.float64(0.503), 'vx': np.float64(0.684), 'vy': np.float64(0.704), 'ax': np.float64(0.483), 'ay': np.float64(0.529)}
static rEFH map, R2 on the test bins: {'x': np.float64(0.189), 'y': np.float64(0.38), 'vx': np.float64(0.594), 'vy': np.float64(0.598), 'ax': np.float64(0.326), 'ay': np.float64(0.37)}
def summary(snr):
# SNR averaged over x and y per kinematic variable (Results: summaries average x and y)
return pd.DataFrame({"pos": snr[["x", "y"]].mean(axis=1), "vel": snr[["vx", "vy"]].mean(axis=1),
"acc": snr[["ax", "ay"]].mean(axis=1)})
S = summary(snr)
print(S.round(2).to_string())
pos vel acc linear 0.80 1.39 0.39 kfobs 5.00 3.30 1.02 kflatent 4.08 3.19 0.98 kflatent_kf 4.26 3.40 1.08 wf 2.96 3.25 1.34 ukf 5.63 3.82 1.11 refh 1.49 3.94 1.86 refh_kf 5.90 5.16 2.37
KFlatent's EM ran all 100 epochs: the last epoch raised the log-likelihood by 4.25e-04 of the total gain, above the Methods' 1/100 of a percent (1e-4), so the 100-epoch limit stopped it, as the Methods allow. The static rEFH map fits its training bins better than the test bins (R2 0.39 and 0.51 against 0.19 and 0.37 for x and y position, 0.67 and 0.71 against 0.59 and 0.59 for velocity); the second KF of the rEFH (+KF) is fitted on these training outputs, as the paper's dynamic map is.
Variants run for every claim of this notebook¶
The rEFH and the EM of KFlatent are stochastic and the paper reports one training run per session; seed 0 is the main run (a standard choice, as the claims file says). Two more rEFH seeds, and two more seeds of the factor-analysis initialization of KFlatent's EM, show how much a claim depends on the run. The code's kinematics (cursor position, no filter, bin means) are the variant the claims file names.
variants = {}
for s in (1, 2):
sp_s = M.prepare(d)
t = time.time()
r = M.dec_refh(sp_s, seed=10 * s + 1)
v = r2.copy()
v.loc["refh"], v.loc["refh_kf"] = sp_s.score(r["refh"]), sp_s.score(r["refh_kf"])
variants[f"rEFH seed {s}"] = M.snr(v)
print(f"rEFH seed {s}: {time.time() - t:.0f} s")
for s in (1, 2):
sp_s = M.prepare(d)
t = time.time()
st, dyn, kfi_s = M.dec_kflatent(sp_s, seed=10 * s, return_info=True)
v = r2.copy()
v.loc["kflatent"], v.loc["kflatent_kf"] = sp_s.score(st), sp_s.score(dyn)
variants[f"KFlatent seed {s}"] = M.snr(v)
print(f"KFlatent seed {s}: {kfi_s['epochs']} E steps, last gain {kfi_s['last_rel_gain']:.2e}, {time.time() - t:.0f} s")
sp_c, pred_c, r2_c, secs, info_c = fit_config(seed=0, reading="code")
variants["code kinematics"] = M.snr(r2_c)
print(f"code kinematics: {secs:.0f} s; KFlatent EM {info_c['kflatent']['epochs']} E steps, last gain {info_c['kflatent']['last_rel_gain']:.2e}")
for k, v in variants.items():
print(f"\n{k}, SNR averaged over x and y:"); print(summary(v).round(2).to_string())
rEFH seed 1: 39 s
rEFH seed 2: 37 s
KFlatent seed 1: 100 E steps, last gain 6.63e-04, 19 s
KFlatent seed 2: 100 E steps, last gain 6.61e-04, 19 s
code kinematics: 58 s; KFlatent EM 100 E steps, last gain 4.25e-04
rEFH seed 1, SNR averaged over x and y:
pos vel acc
linear 0.80 1.39 0.39
kfobs 5.00 3.30 1.02
kflatent 4.08 3.19 0.98
kflatent_kf 4.26 3.40 1.08
wf 2.96 3.25 1.34
ukf 5.63 3.82 1.11
refh 1.53 3.95 1.96
refh_kf 5.93 5.06 2.36
rEFH seed 2, SNR averaged over x and y:
pos vel acc
linear 0.80 1.39 0.39
kfobs 5.00 3.30 1.02
kflatent 4.08 3.19 0.98
kflatent_kf 4.26 3.40 1.08
wf 2.96 3.25 1.34
ukf 5.63 3.82 1.11
refh 1.52 4.07 1.88
refh_kf 5.89 5.20 2.38
KFlatent seed 1, SNR averaged over x and y:
pos vel acc
linear 0.80 1.39 0.39
kfobs 5.00 3.30 1.02
kflatent 4.12 3.18 0.97
kflatent_kf 4.30 3.38 1.08
wf 2.96 3.25 1.34
ukf 5.63 3.82 1.11
refh 1.49 3.94 1.86
refh_kf 5.90 5.16 2.37
KFlatent seed 2, SNR averaged over x and y:
pos vel acc
linear 0.80 1.39 0.39
kfobs 5.00 3.30 1.02
kflatent 4.10 3.22 0.98
kflatent_kf 4.28 3.39 1.09
wf 2.96 3.25 1.34
ukf 5.63 3.82 1.11
refh 1.49 3.94 1.86
refh_kf 5.90 5.16 2.37
code kinematics, SNR averaged over x and y:
pos vel acc
linear 0.73 1.11 0.30
kfobs 4.84 2.95 0.72
kflatent 3.99 3.00 0.79
kflatent_kf 4.28 3.44 1.01
wf 2.88 3.30 1.36
ukf 5.36 3.34 0.78
refh 1.39 3.54 1.64
refh_kf 5.85 5.03 2.22
units-monkey-i¶
After units firing below 0.5 Hz were discarded, monkey I's sessions had between 171 and 413 units, the range reflecting that only 7 of 37 sessions included the S1 array.
rate_session = n_spk / dur
n_cut = int((rate_session > 0.5).sum())
print(f"units in the file {len(ids)}; above 0.5 Hz over the session {n_cut}; at or below 0.5 Hz {len(ids) - n_cut}")
print(f"criterion range 153.9 to 454.3; inside: {153.9 <= n_cut <= 454.3}")
print("hash units in the file: 0 (no unit id ends in 00); channels with sorted units:", len(set(int(i) // 100 for i in ids)))
for thr in (0.25, 0.5, 1.0):
print(f" threshold {thr} Hz: {(rate_session > thr).sum()} units")
units in the file 130; above 0.5 Hz over the session 111; at or below 0.5 Hz 19 criterion range 153.9 to 454.3; inside: False hash units in the file: 0 (no unit id ends in 00); channels with sorted units: 82 threshold 0.25 Hz: 117 units threshold 0.5 Hz: 111 units threshold 1.0 Hz: 94 units
Result: not reproduced. 111 of the 130 units in the file fire above 0.5 Hz over the session, below the 153.9 to 454.3 range. The file holds sorted units only (no unit id ends in 00), on 82 channels. The paper's counts include an unsorted hash unit for each channel (Methods, Data acquisition), up to 96 more here, so the shortfall is consistent with the missing hash units.
fig2-velocity-boost¶
Relative to the static linear map, the rEFH (+KF) raises velocity SNR by more than twice as much as KFobs and nearly twice as much as the UKF.
def velocity_boost(snr):
S = summary(snr)
g = S["vel"] - S.loc["linear", "vel"]
rk, ru = g["refh_kf"] / g["kfobs"], g["refh_kf"] / g["ukf"]
return dict(gain_refh_kf=g["refh_kf"], gain_kfobs=g["kfobs"], gain_ukf=g["ukf"], ratio_kfobs=rk, ratio_ukf=ru,
kfobs_part=rk > 2, ukf_part=1.5 <= ru <= 2.2)
rows = {"main": velocity_boost(snr), **{k: velocity_boost(v) for k, v in variants.items()}}
print(pd.DataFrame(rows).T.round(2).to_string())
gain_refh_kf gain_kfobs gain_ukf ratio_kfobs ratio_ukf kfobs_part ukf_part main 3.762508 1.908801 2.421105 1.971137 1.554046 False True rEFH seed 1 3.668245 1.908801 2.421105 1.921754 1.515112 False True rEFH seed 2 3.804582 1.908801 2.421105 1.99318 1.571424 False True KFlatent seed 1 3.762508 1.908801 2.421105 1.971137 1.554046 False True KFlatent seed 2 3.762508 1.908801 2.421105 1.971137 1.554046 False True code kinematics 3.916383 1.845043 2.229022 2.122651 1.756996 True True
Result: not reproduced, resting on the Methods reading of the kinematics. The rEFH (+KF) raises velocity SNR over linear regression by 3.67 dB, KFobs by 1.91 dB and the UKF by 2.42 dB. That is 1.92 times KFobs's gain, where the criterion needs more than 2, and 1.51 times the UKF's, inside 1.5 to 2.2. With the code's kinematics the ratios are 2.05 and 1.70 and the claim holds. The other two rEFH seeds give KFobs ratios of 1.93 and 1.95.
fig2-kflatent-dynamic-little-effect¶
Adding the dynamic mapping (a second Kalman filter) to KFlatent changes its reconstruction little.
def kfl_dyn(snr):
S = summary(snr)
dd = S.loc["kflatent_kf"] - S.loc["kflatent"]
return {**{f"diff_{k}": v for k, v in dd.items()}, "all_below_0.5": bool((dd.abs() < 0.5).all())}
rows = {"main": kfl_dyn(snr), **{k: kfl_dyn(v) for k, v in variants.items() if k == "code kinematics" or k.startswith("KFlatent")}}
print(pd.DataFrame(rows).T.round(2).to_string())
print("per axis, main (dB):"); print((snr.loc["kflatent_kf"] - snr.loc["kflatent"]).round(2).to_string())
diff_pos diff_vel diff_acc all_below_0.5 main 0.177915 0.20691 0.09447 True KFlatent seed 1 0.175942 0.207203 0.106246 True KFlatent seed 2 0.188465 0.172214 0.103056 True code kinematics 0.28865 0.445339 0.215826 True per axis, main (dB): x 0.17 y 0.18 vx 0.23 vy 0.18 ax 0.08 ay 0.10
Result: reproduced. The dynamic map changes KFlatent's SNR by 0.18 dB for position, 0.21 for velocity and 0.09 for acceleration, all below 0.5 dB. With the other two EM seeds the changes are 0.18 to 0.19, 0.17 to 0.21 and 0.10 to 0.11 dB, and with the code's kinematics 0.29, 0.45 and 0.22 dB.
fig2-refh-static-vel-acc¶
Under the static mapping alone, velocity and acceleration are reconstructed better from the rEFH than from any other filter.
OTHERS = ["linear", "kfobs", "kflatent", "kflatent_kf", "wf", "ukf"]
def static_refh(snr):
S = summary(snr)
out = {}
for v in ("vel", "acc"):
best = S.loc[OTHERS, v].idxmax()
out[f"{v}_refh"] = S.loc["refh", v]; out[f"{v}_best_other"] = f"{best} {S.loc[best, v]:.2f}"
out[f"{v}_margin"] = S.loc["refh", v] - S.loc[best, v]
out["holds"] = bool(out["vel_margin"] > 0 and out["acc_margin"] > 0)
return out
rows = {"main": static_refh(snr), **{k: static_refh(v) for k, v in variants.items()}}
print(pd.DataFrame(rows).T.to_string())
vel_refh vel_best_other vel_margin acc_refh acc_best_other acc_margin holds main 3.940508 ukf 3.82 0.124582 1.862058 wf 1.34 0.526171 True rEFH seed 1 3.949273 ukf 3.82 0.133347 1.959435 wf 1.34 0.623548 True rEFH seed 2 4.073236 ukf 3.82 0.25731 1.880042 wf 1.34 0.544155 True KFlatent seed 1 3.940508 ukf 3.82 0.124582 1.862058 wf 1.34 0.526171 True KFlatent seed 2 3.940508 ukf 3.82 0.124582 1.862058 wf 1.34 0.526171 True code kinematics 3.54417 kflatent_kf 3.44 0.102403 1.63548 wf 1.36 0.271018 True
Result: reproduced. Among the decoders other than the rEFH (+KF), the static rEFH has the highest velocity SNR (3.87 dB against the UKF's 3.82) and acceleration SNR (1.90 against the WF's 1.34). The velocity margin is small: 0.05 dB here, 0.08 and 0.17 dB with the other seeds, and 0.09 dB over KFlatent (+KF) with the code's kinematics. The claim holds in all six runs (main, two rEFH seeds, two EM seeds, code kinematics).
fig2-position-dynamic-gain¶
Position is reconstructed much more accurately from the rEFH with the dynamic than with the static mapping, by about 3 dB of SNR.
def pos_gain(snr):
S = summary(snr)
g = S.loc["refh_kf", "pos"] - S.loc["refh", "pos"]
return dict(static=S.loc["refh", "pos"], dynamic=S.loc["refh_kf", "pos"], gain=g,
gain_x=snr.loc["refh_kf", "x"] - snr.loc["refh", "x"], gain_y=snr.loc["refh_kf", "y"] - snr.loc["refh", "y"],
within_3_pm_0_5=abs(g - 3) <= 0.5)
rows = {"main": pos_gain(snr), **{k: pos_gain(v) for k, v in variants.items()}}
print(pd.DataFrame(rows).T.round(2).to_string())
static dynamic gain gain_x gain_y within_3_pm_0_5 main 1.492066 5.898745 4.406679 4.360466 4.452892 False rEFH seed 1 1.526737 5.926047 4.39931 4.42352 4.375099 False rEFH seed 2 1.517594 5.894416 4.376823 4.315187 4.438458 False KFlatent seed 1 1.492066 5.898745 4.406679 4.360466 4.452892 False KFlatent seed 2 1.492066 5.898745 4.406679 4.360466 4.452892 False code kinematics 1.393767 5.845619 4.451852 4.437322 4.466382 False
Result: not reproduced. The dynamic mapping raises the rEFH's position SNR from 1.46 to 5.91 dB, a gain of 4.45 dB (4.32 for x, 4.58 for y), outside 3 +- 0.5 dB. The gain is larger than the paper's, from 4.22 to 4.52 dB across the other seeds and the code's kinematics. The static map's position R2 is 0.39 and 0.51 on the training bins and 0.19 and 0.37 on the test bins (printed after the fits), so the second KF is fitted on outputs less noisy than those it filters at test.
fig2-sagittal-better¶
For each decoder and each derivative order, the sagittal (y) component is decoded better than the lateral (x) component.
def sagittal(snr):
diff = pd.DataFrame({"pos": snr["y"] - snr["x"], "vel": snr["vy"] - snr["vx"], "acc": snr["ay"] - snr["ax"]})
return diff, int((diff > 0).values.sum())
diff, n = sagittal(snr)
print("SNR y minus x (dB), main:"); print(diff.round(2).to_string())
print(f"pairs with y > x: {n} of 24")
for k, v in variants.items():
dv, nv = sagittal(v)
print(f"{k}: {nv} of 24; pairs with x >= y: {[(i, c) for i in dv.index for c in dv.columns if dv.loc[i, c] <= 0]}")
SNR y minus x (dB), main:
pos vel acc
linear 0.72 0.01 0.29
kfobs 1.12 0.51 0.58
kflatent 1.00 0.33 0.55
kflatent_kf 1.00 0.28 0.57
wf 0.42 0.34 0.54
ukf 1.40 0.82 0.66
refh 1.17 0.04 0.29
refh_kf 1.26 0.42 0.36
pairs with y > x: 24 of 24
rEFH seed 1: 23 of 24; pairs with x >= y: [('refh', 'vel')]
rEFH seed 2: 24 of 24; pairs with x >= y: []
KFlatent seed 1: 24 of 24; pairs with x >= y: []
KFlatent seed 2: 24 of 24; pairs with x >= y: []
code kinematics: 24 of 24; pairs with x >= y: []
Result: partially reproduced. The sagittal component has the higher SNR in 23 of the 24 decoder and variable pairs. The exception is the static rEFH's velocity, where y minus x is -0.06 dB. With rEFH seed 2 all 24 pairs hold; with rEFH seed 1, the two other EM seeds and the code's kinematics 23 hold, with the same exception. Linear regression's velocity margin is 0.01 dB.
fig2-linear-worst¶
Static linear regression reconstructs much worse than the other decoders.
def linear_worst(snr):
m = snr.drop(index="linear") - snr.loc["linear"]
return m, bool((m > 0).values.all())
m, holds = linear_worst(snr)
print("SNR of each decoder minus linear regression (dB), main:"); print(m.round(2).to_string())
print("every decoder above linear regression for every variable:", holds, "; smallest margin", round(m.values.min(), 2))
for k, v in variants.items():
mv, hv = linear_worst(v)
print(f"{k}: holds {hv}, smallest margin {mv.values.min():.2f} dB ({mv.stack().idxmin()})")
SNR of each decoder minus linear regression (dB), main:
x y vx vy ax ay
kfobs 4.00 4.41 1.66 2.16 0.48 0.78
kflatent 3.14 3.42 1.64 1.96 0.46 0.72
kflatent_kf 3.31 3.60 1.87 2.14 0.54 0.83
wf 2.31 2.02 1.69 2.02 0.82 1.07
ukf 4.49 5.17 2.02 2.83 0.54 0.90
refh 0.47 0.92 2.53 2.56 1.47 1.48
refh_kf 4.83 5.37 3.56 3.97 1.95 2.02
every decoder above linear regression for every variable: True ; smallest margin 0.46
rEFH seed 1: holds True, smallest margin 0.46 dB (('kflatent', 'ax'))
rEFH seed 2: holds True, smallest margin 0.46 dB (('kflatent', 'ax'))
KFlatent seed 1: holds True, smallest margin 0.47 dB (('refh', 'x'))
KFlatent seed 2: holds True, smallest margin 0.47 dB (('refh', 'x'))
code kinematics: holds True, smallest margin 0.20 dB (('kfobs', 'ax'))
Result: reproduced. Every decoder has a higher SNR than linear regression for each of the six variables. The smallest margin is 0.45 dB (static rEFH, x position), and the claim holds in every variant (smallest margin 0.20 dB, with the code's kinematics).
fig2-gross-failures¶
Gross failures (negative coefficients of determination) occur only for the UKF and KFlatent (+KF), in a few sessions.
def gross(r2):
neg = [(i, c, round(r2.loc[i, c], 3)) for i in r2.index for c in r2.columns if r2.loc[i, c] < 0]
return neg, all(i in ("ukf", "kflatent_kf") for i, _, _ in neg)
neg, holds = gross(r2)
print("negative R2, main:", neg, "; only UKF or KFlatent (+KF):", holds)
print("lowest R2 per decoder, main:"); print(r2.min(axis=1).round(3).to_string())
neg_c, holds_c = gross(r2_c)
print("negative R2, code kinematics:", neg_c, holds_c)
negative R2, main: [] ; only UKF or KFlatent (+KF): True lowest R2 per decoder, main: linear 0.055 kfobs 0.154 kflatent 0.150 kflatent_kf 0.166 wf 0.217 ukf 0.165 refh 0.189 refh_kf 0.397 negative R2, code kinematics: [] True
Result: reproduced, with no failures to count. No decoder has a negative R2 at 64 ms in this session; the lowest is linear regression's 0.055. The criterion, that every negative R2 belongs to the UKF or KFlatent (+KF), therefore holds trivially. Whether failures occur in a few sessions cannot be tested on one session.
fig4-vel-acc-every-session¶
For velocity and acceleration, the rEFH (+KF) is better than KFobs and the UKF in every one of the 196 session, variable and axis combinations, and significantly so across sessions.
def every(snr):
cols = ["vx", "vy", "ax", "ay"]
m = pd.DataFrame({r: snr.loc["refh_kf", cols] - snr.loc[r, cols] for r in ("kfobs", "ukf")}).T
return m, int((m > 0).values.sum())
m, n = every(snr)
print("SNR rEFH (+KF) minus rival (dB), main:"); print(m.round(2).to_string())
print(f"comparisons favouring the rEFH (+KF): {n} of 8")
for k, v in variants.items():
mv, nv = every(v)
print(f"{k}: {nv} of 8, smallest margin {mv.values.min():.2f} dB")
SNR rEFH (+KF) minus rival (dB), main:
vx vy ax ay
kfobs 1.90 1.81 1.47 1.24
ukf 1.54 1.14 1.41 1.11
comparisons favouring the rEFH (+KF): 8 of 8
rEFH seed 1: 8 of 8, smallest margin 1.08 dB
rEFH seed 2: 8 of 8, smallest margin 1.19 dB
KFlatent seed 1: 8 of 8, smallest margin 1.11 dB
KFlatent seed 2: 8 of 8, smallest margin 1.11 dB
code kinematics: 8 of 8, smallest margin 1.28 dB
Result: reproduced for this session. The rEFH (+KF) beats KFobs and the UKF in all 8 velocity and acceleration comparisons, by 1.04 to 1.81 dB, and in every variant (smallest margin 1.04 dB). The bootstrap across the paper's 49 sessions cannot be run on one session.
fig4-best-variable¶
The best-decoded kinematic variable is velocity for the rEFH (+KF), but position for KFobs and the UKF.
def best_var(snr):
S = summary(snr)
b = S.loc[["refh_kf", "kfobs", "ukf"]].idxmax(axis=1)
return S.loc[["refh_kf", "kfobs", "ukf"]], b, bool(b["refh_kf"] == "vel" and b["kfobs"] == "pos" and b["ukf"] == "pos")
Sb, b, holds = best_var(snr)
print(Sb.round(2).to_string()); print("best variable:", b.to_dict(), "; holds:", holds)
for k, v in variants.items():
_, bv, hv = best_var(v)
print(f"{k}: {bv.to_dict()} holds {hv}")
pos vel acc
refh_kf 5.90 5.16 2.37
kfobs 5.00 3.30 1.02
ukf 5.63 3.82 1.11
best variable: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} ; holds: False
rEFH seed 1: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} holds False
rEFH seed 2: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} holds False
KFlatent seed 1: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} holds False
KFlatent seed 2: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} holds False
code kinematics: {'refh_kf': 'pos', 'kfobs': 'pos', 'ukf': 'pos'} holds False
Result: not reproduced. Position is the best-decoded variable for all three decoders. The rEFH (+KF) has 5.91 dB for position against 5.06 for velocity; KFobs and the UKF have position highest, as the paper says. The ordering is the same with the other seeds and the code's kinematics.
Single-session parts of claims blocked for the other sessions¶
The claims below need a bootstrap or a regression over the paper's 49 sessions, which the dandiset does not hold, so they stay blocked. Their point comparisons can be computed on this one session and are printed here; they do not decide an outcome.
fig2-refh-kf-best¶
Averaged over the 49 sessions, no decoder outperforms the rEFH (+KF) for any kinematic variable, and for velocity and acceleration it is clearly superior to all other decoders.
def refh_best(snr):
S = summary(snr)
m = S.loc["refh_kf"] - S.drop(index="refh_kf").max()
return {**{f"{v}_refh_kf": S.loc["refh_kf", v] for v in S.columns},
**{f"{v}_margin_over_next": m[v] for v in S.columns},
**{f"{v}_next": S.drop(index="refh_kf")[v].idxmax() for v in S.columns}}
rows = {"main": refh_best(snr), **{k: refh_best(v) for k, v in variants.items()}}
print(pd.DataFrame(rows).T.to_string())
pos_refh_kf vel_refh_kf acc_refh_kf pos_margin_over_next vel_margin_over_next acc_margin_over_next pos_next vel_next acc_next main 5.898745 5.157329 2.373325 0.27042 1.216822 0.511267 ukf refh refh rEFH seed 1 5.926047 5.063067 2.357734 0.297722 1.113794 0.398299 ukf refh refh rEFH seed 2 5.894416 5.199404 2.381027 0.266092 1.126168 0.500985 ukf refh refh KFlatent seed 1 5.898745 5.157329 2.373325 0.27042 1.216822 0.511267 ukf refh refh KFlatent seed 2 5.898745 5.157329 2.373325 0.27042 1.216822 0.511267 ukf refh refh code kinematics 5.845619 5.025021 2.215099 0.489548 1.480851 0.57962 ukf refh refh
intro-gain-over-nearest¶
The rEFH (+KF) improves SNR over its nearest competitor on average by 0.14 dB for position, 1.60 dB for velocity and 1.2 dB for acceleration; the position gain is not significant at p = 0.01.
def gains(snr):
S = summary(snr)
a = S.loc["refh_kf"] - S.loc[OTHERS].max()
b = S.loc["refh_kf"] - S.loc["ukf"]
return {**{f"a_{v}": a[v] for v in S.columns}, **{f"b_{v}": b[v] for v in S.columns},
**{f"a_rival_{v}": S.loc[OTHERS, v].idxmax() for v in S.columns}}
rows = {"main": gains(snr), **{k: gains(v) for k, v in variants.items()}}
print("expected 0.14 +- 0.014, 1.6 +- 0.16, 1.2 +- 0.12 dB (pos, vel, acc)")
print(pd.DataFrame(rows).T.to_string())
expected 0.14 +- 0.014, 1.6 +- 0.16, 1.2 +- 0.12 dB (pos, vel, acc)
a_pos a_vel a_acc b_pos b_vel b_acc a_rival_pos a_rival_vel a_rival_acc
main 0.27042 1.341403 1.037437 0.27042 1.341403 1.262569 ukf ukf wf
rEFH seed 1 0.297722 1.24714 1.021847 0.297722 1.24714 1.246979 ukf ukf wf
rEFH seed 2 0.266092 1.383478 1.04514 0.266092 1.383478 1.270272 ukf ukf wf
KFlatent seed 1 0.27042 1.341403 1.037437 0.27042 1.341403 1.262569 ukf ukf wf
KFlatent seed 2 0.27042 1.341403 1.037437 0.27042 1.341403 1.262569 ukf ukf wf
code kinematics 0.489548 1.583254 0.850638 0.489548 1.687361 1.436846 ukf kflatent_kf wf
fig4-position¶
The rEFH (+KF) decodes position better than KFobs, but not distinguishably from the UKF.
rows = {k: dict(minus_kfobs=summary(v).loc["refh_kf", "pos"] - summary(v).loc["kfobs", "pos"],
minus_ukf=summary(v).loc["refh_kf", "pos"] - summary(v).loc["ukf", "pos"])
for k, v in {"main": snr, **variants}.items()}
print("position SNR, rEFH (+KF) minus rival (dB):"); print(pd.DataFrame(rows).T.round(2).to_string())
position SNR, rEFH (+KF) minus rival (dB):
minus_kfobs minus_ukf
main 0.89 0.27
rEFH seed 1 0.92 0.30
rEFH seed 2 0.89 0.27
KFlatent seed 1 0.89 0.27
KFlatent seed 2 0.89 0.27
code kinematics 1.00 0.49
static-refh-position-worse¶
Under the static mapping alone the rEFH reconstructs position significantly worse than the UKF, while the rEFH (+KF) is not distinguishable from it.
rows = {k: dict(static_refh=summary(v).loc["refh", "pos"], ukf=summary(v).loc["ukf", "pos"],
static_refh_minus_ukf=summary(v).loc["refh", "pos"] - summary(v).loc["ukf", "pos"])
for k, v in {"main": snr, **variants}.items()}
print("position SNR (dB):"); print(pd.DataFrame(rows).T.round(2).to_string())
position SNR (dB):
static_refh ukf static_refh_minus_ukf
main 1.49 5.63 -4.14
rEFH seed 1 1.53 5.63 -4.10
rEFH seed 2 1.52 5.63 -4.11
KFlatent seed 1 1.49 5.63 -4.14
KFlatent seed 2 1.49 5.63 -4.14
code kinematics 1.39 5.36 -3.96
Result: blocked, single-session values only. These claims need the bootstrap over the paper's 49 sessions. On this session:
- fig2-refh-kf-best: the rEFH (+KF) has the highest SNR for position, velocity and acceleration (5.91, 5.06 and 2.36 dB), ahead of the next decoder by 0.28 dB (UKF), 1.19 and 0.45 dB (static rEFH); it is highest in every variant.
- intro-gain-over-nearest: the gains over the nearest rival are 0.28, 1.24 and 1.02 dB under reading a (UKF, UKF, WF) and 0.28, 1.24 and 1.25 dB under reading b (UKF). The position and velocity gains are outside 0.14 +- 0.014 and 1.6 +- 0.16 dB under both readings and with every seed (position 0.11 to 0.28 dB, velocity 1.24 to 1.31 dB). With the code's kinematics the velocity gain is inside (1.45 and 1.55 dB) but the position gain is 0.51 dB.
- fig4-position: the rEFH (+KF) position SNR exceeds KFobs's by 0.90 dB and the UKF's by 0.28 dB (0.73 to 1.02 and 0.11 to 0.51 dB across the variants).
- static-refh-position-worse: the static rEFH's position SNR is 4.17 dB below the UKF's (1.46 against 5.63), and 4.01 to 4.12 dB below in the variants.
kflatent-modest¶
KFlatent is only modestly better than the standard Kalman filter; nonlinear tuning (the UKF) improves on the standard filter more.
def kfl_modest(snr):
S = summary(snr)
reads = {"a velocity": ["vel"], "b velocity and acceleration": ["vel", "acc"], "c all three": ["pos", "vel", "acc"]}
out = {}
for name, cols in reads.items():
kfl = max(S.loc["kflatent", cols].mean(), S.loc["kflatent_kf", cols].mean())
kfo, ukf = S.loc["kfobs", cols].mean(), S.loc["ukf", cols].mean()
out[name] = dict(kflatent_better_minus_kfobs=kfl - kfo, ukf_minus_kfobs=ukf - kfo,
holds=bool(kfl > kfo and (ukf - kfo) > (kfl - kfo)))
return pd.DataFrame(out).T
print("main:"); print(kfl_modest(snr).to_string())
for k, v in variants.items():
if k == "code kinematics" or k.startswith("KFlatent"):
print(f"{k}:"); print(kfl_modest(v).to_string())
print("per variable, main (dB):"); print(summary(snr).loc[["kfobs", "kflatent", "kflatent_kf", "ukf"]].round(2).to_string())
main:
kflatent_better_minus_kfobs ukf_minus_kfobs holds
a velocity 0.095547 0.512304 True
b velocity and acceleration 0.077032 0.302999 True
c all three -0.198198 0.409967 False
KFlatent seed 1:
kflatent_better_minus_kfobs ukf_minus_kfobs holds
a velocity 0.08106 0.512304 True
b velocity and acceleration 0.072013 0.302999 True
c all three -0.187437 0.409967 False
KFlatent seed 2:
kflatent_better_minus_kfobs ukf_minus_kfobs holds
a velocity 0.084971 0.512304 True
b velocity and acceleration 0.077652 0.302999 True
c all three -0.188473 0.409967 False
code kinematics:
kflatent_better_minus_kfobs ukf_minus_kfobs holds
a velocity 0.488086 0.383979 False
b velocity and acceleration 0.389067 0.222769 False
c all three 0.071948 0.319981 True
per variable, main (dB):
pos vel acc
kfobs 5.00 3.30 1.02
kflatent 4.08 3.19 0.98
kflatent_kf 4.26 3.40 1.08
ukf 5.63 3.82 1.11
Result: reproduced. On velocity (reading a), the better KFlatent decoder exceeds KFobs by 0.10 dB while the UKF exceeds it by 0.51 dB; on velocity and acceleration (reading b) by 0.08 and 0.30 dB. Over all three variables (reading c) KFlatent is 0.20 dB below KFobs. With the two other EM seeds KFlatent's lead over KFobs is 0.08 and 0.08 dB under reading a and 0.07 and 0.08 dB under reading b, so readings a and b hold for all three EM runs. With the code's kinematics readings a and b fail (KFlatent ahead of KFobs by 0.49 and 0.39 dB, more than the UKF's 0.38 and 0.22) and reading c holds (0.07 against 0.32 dB). Under either kinematics the claim holds under one reading.
fig8-velocity-variants¶
Under the static mapping, about 0.7 dB of the rEFH's velocity advantage over the UKF is lost when the updated visible units are left out of decoding, and 0.4 dB when standard-normal visible units replace Poisson ones; Kalman filtering does not remove these differences.
The variant without visible units decodes the canonical rEFH from its hidden means alone. The standard-normal rEFH is trained on standardized counts (the Results, "Variant decoders"); the code's transform 2 sqrt(R + 3/8) is the variant the criterion names.
t = time.time()
sp8 = M.prepare(d)
canon = M.dec_refh(sp8, seed=1, variants=True) # seed 1 is the main rEFH run (run_all uses seed + 1)
norm_std = M.dec_refh(sp8, seed=1, units="normal", variants="standardize")
norm_ans = M.dec_refh(sp8, seed=1, units="normal", variants="anscombe")
print(f"variants fitted in {time.time() - t:.0f} s")
v8 = pd.DataFrame({
"canonical": sp8.score(canon["refh"]), "canonical +KF": sp8.score(canon["refh_kf"]),
"hidden only": sp8.score(canon["refh_hid"]), "hidden only +KF": sp8.score(canon["refh_hid_kf"]),
"std normal (standardized)": sp8.score(norm_std["refh"]), "std normal (standardized) +KF": sp8.score(norm_std["refh_kf"]),
"std normal (code transform)": sp8.score(norm_ans["refh"]), "std normal (code transform) +KF": sp8.score(norm_ans["refh_kf"]),
"ukf": r2.loc["ukf"].values}, index=M.VARS).T
S8 = summary(M.snr(v8))
print(S8.round(2).to_string())
vel = S8["vel"]
loss_hid = vel["canonical"] - vel["hidden only"]
loss_std = vel["canonical"] - vel["std normal (standardized)"]
loss_ans = vel["canonical"] - vel["std normal (code transform)"]
print(f"\nstatic velocity loss, hidden only: {loss_hid:.2f} dB (text 0.7 +- 0.07; figure about 0.3)")
print(f"static velocity loss, standard normal (standardized): {loss_std:.2f} dB (text 0.4 +- 0.04; figure about 0.4)")
print(f"static velocity loss, standard normal (code transform): {loss_ans:.2f} dB")
print("advantage of canonical static over UKF, velocity:", round(vel["canonical"] - vel["ukf"], 2))
for nm in ["hidden only +KF", "std normal (standardized) +KF", "std normal (code transform) +KF"]:
print(f"dynamic velocity, canonical +KF minus {nm}: {vel['canonical +KF'] - vel[nm]:.2f} dB")
variants fitted in 79 s
pos vel acc
canonical 1.49 3.94 1.86
canonical +KF 5.90 5.16 2.37
hidden only 1.43 3.73 1.65
hidden only +KF 5.75 4.91 2.21
std normal (standardized) 1.41 3.89 1.84
std normal (standardized) +KF 5.82 4.89 2.24
std normal (code transform) 0.99 1.71 0.67
std normal (code transform) +KF 4.75 3.44 1.24
ukf 5.63 3.82 1.11
static velocity loss, hidden only: 0.21 dB (text 0.7 +- 0.07; figure about 0.3)
static velocity loss, standard normal (standardized): 0.05 dB (text 0.4 +- 0.04; figure about 0.4)
static velocity loss, standard normal (code transform): 2.23 dB
advantage of canonical static over UKF, velocity: 0.12
dynamic velocity, canonical +KF minus hidden only +KF: 0.25 dB
dynamic velocity, canonical +KF minus std normal (standardized) +KF: 0.27 dB
dynamic velocity, canonical +KF minus std normal (code transform) +KF: 1.72 dB
Result: not reproduced. Leaving the updated visible units out costs 0.11 dB of static velocity SNR (the text says 0.7 dB, figure 8 about 0.3). The standard-normal rEFH trained on standardized counts gains 0.08 dB (a loss of -0.08 dB against the text's 0.4). Trained on the code's transform, 2 sqrt(R + 3/8), it loses 2.11 dB. Under the dynamic mapping the canonical rEFH (+KF) still exceeds the variants, by 0.15, 0.15 and 1.68 dB. The static rEFH's own velocity advantage over the UKF is only 0.05 dB here.
fig8-gaussian-still-beats-ukf¶
The rEFH keeps an advantage over the UKF when its visible units are standard normal instead of Poisson, so the noise model accounts for only part of its advantage.
for nm in ["std normal (standardized) +KF", "std normal (code transform) +KF"]:
dv = S8.loc[nm] - S8.loc["ukf"]
print(f"{nm} minus UKF (dB): vel {dv['vel']:.2f}, acc {dv['acc']:.2f}, pos {dv['pos']:.2f}; vel and acc above UKF: {bool(dv['vel'] > 0 and dv['acc'] > 0)}")
std normal (standardized) +KF minus UKF (dB): vel 1.07, acc 1.13, pos 0.19; vel and acc above UKF: True std normal (code transform) +KF minus UKF (dB): vel -0.37, acc 0.13, pos -0.88; vel and acc above UKF: False
Result: reproduced, resting on the standardized-counts reading. Trained on standardized counts, the standard-normal rEFH (+KF) exceeds the UKF by 1.09 dB for velocity and 1.17 dB for acceleration. Trained on the code's transform it is 0.43 dB below the UKF for velocity (0.12 dB above for acceleration), and the claim fails.
fig8-position-variants (single-session part)¶
For position the rEFH variants change the static but not the dynamic mapping; with the Kalman filter all three variants decode position about equally.
The dynamic part needs the bootstrap over sessions; the static part and the dynamic differences are printed.
pos = S8["pos"]
print(pos.round(2).to_string())
print(f"static: canonical minus hidden only {pos['canonical'] - pos['hidden only']:.2f} dB, minus standard normal (standardized) "
f"{pos['canonical'] - pos['std normal (standardized)']:.2f} dB, minus standard normal (code transform) {pos['canonical'] - pos['std normal (code transform)']:.2f} dB")
for nm in ["hidden only +KF", "std normal (standardized) +KF", "std normal (code transform) +KF"]:
print(f"dynamic: canonical +KF minus {nm}: {pos['canonical +KF'] - pos[nm]:.2f} dB")
canonical 1.49 canonical +KF 5.90 hidden only 1.43 hidden only +KF 5.75 std normal (standardized) 1.41 std normal (standardized) +KF 5.82 std normal (code transform) 0.99 std normal (code transform) +KF 4.75 ukf 5.63 static: canonical minus hidden only 0.06 dB, minus standard normal (standardized) 0.09 dB, minus standard normal (code transform) 0.51 dB dynamic: canonical +KF minus hidden only +KF: 0.15 dB dynamic: canonical +KF minus std normal (standardized) +KF: 0.08 dB dynamic: canonical +KF minus std normal (code transform) +KF: 1.15 dB
Result: blocked, single-session values only. Under the static mapping the canonical rEFH has a higher position SNR (1.46 dB) than both variants (hidden units only 1.43, standard normal on standardized counts 1.38, on the code's transform 0.99), so the static part holds. Under the dynamic mapping the canonical rEFH (+KF) is 0.11 and 0.16 dB above the two variants (1.17 dB above the code-transform variant); whether that is below p < 0.01 needs the bootstrap over sessions.
Summary¶
| Claim | Outcome |
|---|---|
| units-monkey-i | not reproduced (111 units; no hash units in the file) |
| fig2-velocity-boost | not reproduced, resting on the Methods kinematics (reproduced with the code's) |
| fig2-kflatent-dynamic-little-effect | reproduced |
| fig2-refh-static-vel-acc | reproduced |
| fig2-position-dynamic-gain | not reproduced (gain 4.45 dB, not about 3) |
| fig2-sagittal-better | partially reproduced (23 of 24) |
| fig2-linear-worst | reproduced |
| fig2-gross-failures | reproduced (no negative R2 at 64 ms; other configurations have some, see binwidth_training_neurons.ipynb) |
| fig4-vel-acc-every-session | reproduced for this session (bootstrap not run) |
| fig4-best-variable | not reproduced (the rEFH (+KF) decodes position best) |
| kflatent-modest | reproduced |
| fig8-velocity-variants | not reproduced |
| fig8-gaussian-still-beats-ukf | reproduced, resting on standardized counts |
| fig5-advantage-with-binwidth, fig5-derivative-order, advantage-grows-with-derivative, fig7-affine-position, fig7-kinematic-model-saturates | see binwidth_training_neurons.ipynb |
| fig2-refh-kf-best, intro-gain-over-nearest, fig4-advantage-grows, fig4-position, fig5-vel-acc-all-bins, fig5-position-bins, fig6-training-time, fig7-vel-acc-neurons, fig7-position-neurons, fig8-position-variants, fig9-taps, static-refh-position-worse | blocked: they need the paper's other sessions, which the dandiset does not hold (single-session values printed above and in binwidth_training_neurons.ipynb) |
| fig4-task-similar | blocked: monkey C's sessions have no archive identifier, and the dandiset's one monkey I session is probably not one of the paper's |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. One session in place of 49, so every mean over sessions is this session's value and no bootstrap is run. The session is dated 2017-02-02, which is not in the paper's session list. Sorted units only (111 above 0.5 Hz), without the hash units the paper kept. The test part is the last 329 s of the training file; the test file's kinematics are withheld. Six bins that overlap the file's three 0.1 s gaps are left out of fitting and scoring. All decoders are re-implemented in Python from the Methods and the paper's code; the UKF re-implements Li's UKF2, whose fitting code is not in the paper's repository. The rEFH is trained in single precision. EM for KFlatent ran the full 100 epochs without meeting the convergence criterion of 1/100 of a percent, as the Methods allow. The rEFH and the factor-analysis initialization of EM are each run with three seeds. The Wiener filter takes the history of the first test bins from the last training bins.
What was found about the dataset. The dandiset holds one of the paper's three monkeys and one session, which the paper's session list does not include. It has no hash units. Three 0.1 s gaps split the training file into four segments, with NaN kinematics and no spikes, and nothing says what was cut. The test file's 190 s have no kinematics. The units and the kinematics share one clock: velocity is best decoded from spikes 64 ms earlier.
print(f"notebook run time {(time.time() - T0) / 60:.1f} min")
notebook run time 5.3 min