Reproducing Makin et al. (2018) from dandiset 000129: bin width, training data and neurons¶
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, streamed. It is the only published version; the draft holds the same two assets
and differs only in its title. The claims file, claims/000129/10.1088-1741-2552-aa9e95.yaml, was
written from the version of record and the paper's code (jgmakin/rbmish at commit 700062f5).
The dataset holds one session of monkey I (M1 only, 649.1 s with kinematics); the paper used 49
sessions. This notebook refits all eight decoders on that session at bin widths of 16, 32, 64 and
128 ms (figure 5), with 80, 160 and 320 s of training data (figure 6), and with fractions 1/8,
1/4, 1/2 and 1 of the units (figure 7). The data checks and the decoders at 64 ms are in
reproduce_makin_2018.ipynb.
| Claim | Here |
|---|---|
| fig5-advantage-with-binwidth, fig5-derivative-order | tested on one session |
| advantage-grows-with-derivative | tested on one session |
| fig7-affine-position, fig7-kinematic-model-saturates | tested on one session |
| fig5-vel-acc-all-bins, fig5-position-bins, fig6-training-time, fig7-vel-acc-neurons, fig7-position-neurons | blocked: bootstrap across the paper's 49 sessions, which the dandiset does not hold; single-session values printed |
| all others | see reproduce_makin_2018.ipynb |
import os
os.environ.setdefault("OMP_NUM_THREADS", "4")
import time
import numpy as np
import pandas as pd
import makin as M
pd.set_option("display.width", 200)
T0 = time.time()
d = M.load_session()
print(f"dandiset {M.DANDISET} version {M.VERSION}; {d['session_description'][:80]}")
dandiset 000129 version 0.241017.1444; Data from monkey Indy performing self-paced random target reaching task. This fi
Data check¶
print(f"units {len(d['unit_ids'])}; kinematic samples {len(d['finger_pos'])} at {d['rate']:.0f} Hz; NaN samples {np.isnan(d['finger_pos'][:, 0]).sum()}")
rows = []
for bms in (16, 32, 64, 128):
R, X, _ = M.bin_data(d, bms)
ntr = int(320 / (bms / 1000))
keep = M.unit_filter(R, bms, ntr)
rows.append(dict(bin_ms=bms, bins=len(X), train_bins=ntr, test_bins=len(X) - ntr, nan_bins=int(np.isnan(X).any(1).sum()),
units_used=int(keep.sum()), mean_count=R[:, keep].mean().round(3),
speed_max=np.nanmax(np.hypot(X[:, 2], X[:, 3])).round(0)))
print(pd.DataFrame(rows).to_string(index=False))
units 130; kinematic samples 649100 at 1000 Hz; NaN samples 300
bin_ms bins train_bins test_bins nan_bins units_used mean_count speed_max
16 40568 20000 20568 21 111 0.075 546.0
32 20284 10000 10284 12 111 0.150 546.0
64 10142 5000 5142 6 111 0.299 546.0
128 5071 2500 2571 3 111 0.598 512.0
The number of training bins is 320 s divided by the bin width (20,000, 10,000, 5,000 and 2,500, as the Methods give), and the same units pass the 0.5 Hz filter at every width.
Fits¶
Each configuration fits the eight decoders with seed 0 (a standard choice: one training run, as the paper reports one per session). The rEFH is trained for max(floor(Ntrain/250), 20) epochs (80 at 16 ms, 40 at 32 ms, 20 at 64 and 128 ms; Methods and code). For 80 and 160 s of training the paper does not say which bins are tested; the code (getLatentsBMI.m, getTestBins) tests on every bin after the training bins, and this is the main reading; scoring only the bins after 320 s, the same test set as for 320 s, is shown as a variant.
def summary(snr):
return pd.DataFrame({"pos": snr[["x", "y"]].mean(axis=1), "vel": snr[["vx", "vy"]].mean(axis=1),
"acc": snr[["ax", "ay"]].mean(axis=1)})
results, results_320, infos = {}, {}, {}
configs = [("bin", 16, dict(bin_ms=16)), ("bin", 32, dict(bin_ms=32)), ("bin", 64, dict(bin_ms=64)),
("bin", 128, dict(bin_ms=128)), ("train", 80, dict(train_s=80.0)), ("train", 160, dict(train_s=160.0)),
("frac", 1 / 8, dict(fraction=1 / 8)), ("frac", 1 / 4, dict(fraction=1 / 4)), ("frac", 1 / 2, dict(fraction=1 / 2))]
for kind, val, kw in configs:
t = time.time()
sp = M.prepare(d, **kw)
info = {}
pred = M.run_all(sp, seed=0, info=info)
r2 = pd.DataFrame({k: sp.score(v) for k, v in pred.items()}, index=M.VARS).T.loc[M.DECODERS]
results[(kind, val)] = r2
infos[(kind, val)] = info
if kind == "train":
results_320[(kind, val)] = pd.DataFrame({k: sp.score(v, from_s=320.0) for k, v in pred.items()}, index=M.VARS).T.loc[M.DECODERS]
print(f"{kind} {val}: units {sp.Rtr.shape[1]}, train bins {len(sp.Rtr)}, test bins {sp.okte.sum()}, "
f"negative R2 {int((r2.values < 0).sum())}, KFlatent EM {info['kflatent']['epochs']} E steps "
f"(last gain {info['kflatent']['last_rel_gain']:.1e}), {time.time() - t:.0f} s")
results[("train", 320)] = results[("bin", 64)]
results_320[("train", 320)] = results[("bin", 64)]
results[("frac", 1.0)] = results[("bin", 64)]
SNR = {k: summary(M.snr(v)) for k, v in results.items()}
for k, v in SNR.items():
print(f"\n{k}, SNR (dB) averaged over x and y:"); print(v.round(2).to_string())
bin 16: units 111, train bins 20000, test bins 20554, negative R2 2, KFlatent EM 100 E steps (last gain 6.3e-04), 240 s
bin 32: units 111, train bins 10000, test bins 10276, negative R2 0, KFlatent EM 100 E steps (last gain 8.8e-04), 114 s
bin 64: units 111, train bins 5000, test bins 5138, negative R2 0, KFlatent EM 100 E steps (last gain 4.3e-04), 57 s
bin 128: units 111, train bins 2500, test bins 2569, negative R2 0, KFlatent EM 100 E steps (last gain 4.9e-04), 48 s
train 80: units 111, train bins 1250, test bins 8886, negative R2 1, KFlatent EM 100 E steps (last gain 7.2e-04), 44 s
train 160: units 111, train bins 2500, test bins 7636, negative R2 1, KFlatent EM 100 E steps (last gain 4.4e-04), 49 s
frac 0.125: units 13, train bins 5000, test bins 5138, negative R2 13, KFlatent EM 100 E steps (last gain 9.4e-04), 11 s
frac 0.25: units 27, train bins 5000, test bins 5138, negative R2 6, KFlatent EM 100 E steps (last gain 5.3e-04), 13 s
frac 0.5: units 55, train bins 5000, test bins 5138, negative R2 0, KFlatent EM 100 E steps (last gain 6.3e-04), 26 s
('bin', 16), SNR (dB) averaged over x and y:
pos vel acc
linear 0.17 0.46 0.14
kfobs 4.13 2.52 0.46
kflatent 1.88 2.19 0.45
kflatent_kf 3.30 2.60 0.13
wf 3.13 3.30 1.39
ukf 4.68 2.82 0.50
refh 0.18 0.54 0.20
refh_kf 4.42 2.96 0.95
('bin', 32), SNR (dB) averaged over x and y:
pos vel acc
linear 0.42 0.83 0.25
kfobs 4.70 3.01 0.79
kflatent 3.10 2.77 0.67
kflatent_kf 3.80 3.06 0.68
wf 3.00 3.28 1.37
ukf 5.29 3.51 0.92
refh 0.51 1.68 0.64
refh_kf 5.24 4.02 1.69
('bin', 64), SNR (dB) 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.49 3.94 1.86
refh_kf 5.90 5.16 2.37
('bin', 128), SNR (dB) averaged over x and y:
pos vel acc
linear 1.33 2.22 0.51
kfobs 5.19 3.47 1.03
kflatent 4.63 3.43 1.09
kflatent_kf 4.86 3.58 1.22
wf 2.91 3.17 1.24
ukf 5.52 3.86 1.12
refh 3.55 5.00 2.31
refh_kf 6.39 5.45 2.44
('train', 80), SNR (dB) averaged over x and y:
pos vel acc
linear 0.19 0.87 -0.03
kfobs 2.73 2.56 0.60
kflatent 2.02 2.35 0.55
kflatent_kf 3.43 3.06 0.91
wf 1.08 2.13 0.77
ukf 3.41 2.78 0.72
refh 0.44 2.36 0.97
refh_kf 2.48 3.31 1.59
('train', 160), SNR (dB) averaged over x and y:
pos vel acc
linear 0.65 1.20 0.21
kfobs 3.91 2.97 0.85
kflatent 3.32 2.77 0.77
kflatent_kf 3.69 3.11 1.02
wf 2.44 2.83 1.12
ukf 4.27 3.28 0.95
refh 1.07 3.29 1.43
refh_kf 4.36 4.31 1.95
('frac', 0.125), SNR (dB) averaged over x and y:
pos vel acc
linear -0.09 0.27 0.08
kfobs 0.61 0.77 0.18
kflatent -0.10 0.41 0.03
kflatent_kf -2.61 -0.07 -0.21
wf 0.22 0.94 0.29
ukf 0.56 0.71 0.13
refh -0.08 0.31 0.14
refh_kf 0.68 0.80 0.26
('frac', 0.25), SNR (dB) averaged over x and y:
pos vel acc
linear -0.03 0.41 0.17
kfobs 0.32 1.30 0.40
kflatent 0.10 1.10 0.26
kflatent_kf 0.07 1.03 0.15
wf 0.32 1.56 0.53
ukf 0.33 1.31 0.40
refh -0.06 0.63 0.34
refh_kf 0.55 1.56 0.74
('frac', 0.5), SNR (dB) averaged over x and y:
pos vel acc
linear 0.33 0.91 0.24
kfobs 3.11 2.21 0.65
kflatent 1.77 1.93 0.41
kflatent_kf 1.12 1.87 0.40
wf 1.65 2.24 0.65
ukf 3.15 2.34 0.65
refh 0.42 2.01 0.67
refh_kf 3.33 2.83 1.27
('train', 320), SNR (dB) 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.49 3.94 1.86
refh_kf 5.90 5.16 2.37
('frac', 1.0), SNR (dB) 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.49 3.94 1.86
refh_kf 5.90 5.16 2.37
for k, v in results.items():
if k in (("train", 320), ("frac", 1.0)):
continue # the same fit as ("bin", 64)
print(f"\n{k}, R2 per decoder and variable:"); print(v.round(3).to_string())
print("static rEFH map, training R2:", dict(zip(M.VARS, infos[k]["refh_train_r2"].round(3))))
('bin', 16), R2 per decoder and variable:
x y vx vy ax ay
linear -0.000 0.076 0.103 0.098 0.020 0.043
kfobs 0.517 0.690 0.381 0.493 0.009 0.185
kflatent 0.182 0.484 0.362 0.430 0.057 0.138
kflatent_kf 0.433 0.614 0.422 0.478 -0.064 0.116
wf 0.479 0.547 0.513 0.550 0.228 0.317
ukf 0.553 0.740 0.381 0.559 0.016 0.193
refh 0.006 0.073 0.133 0.102 0.036 0.055
refh_kf 0.550 0.710 0.437 0.546 0.116 0.269
static rEFH map, training R2: {'x': np.float64(0.107), 'y': np.float64(0.146), 'vx': np.float64(0.156), 'vy': np.float64(0.156), 'ax': np.float64(0.094), 'ay': np.float64(0.121)}
('bin', 32), R2 per decoder and variable:
x y vx vy ax ay
linear 0.037 0.143 0.174 0.172 0.036 0.077
kfobs 0.597 0.714 0.459 0.538 0.096 0.232
kflatent 0.405 0.598 0.467 0.477 0.103 0.182
kflatent_kf 0.512 0.644 0.493 0.519 0.086 0.199
wf 0.469 0.527 0.511 0.548 0.225 0.315
ukf 0.633 0.762 0.490 0.611 0.116 0.258
refh 0.062 0.158 0.327 0.314 0.111 0.164
refh_kf 0.638 0.753 0.571 0.634 0.265 0.374
static rEFH map, training R2: {'x': np.float64(0.204), 'y': np.float64(0.275), 'vx': np.float64(0.377), 'vy': np.float64(0.385), 'ax': np.float64(0.227), 'ay': np.float64(0.256)}
('bin', 64), R2 per decoder and variable:
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
static rEFH map, training R2: {'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)}
('bin', 128), R2 per decoder and variable:
x y vx vy ax ay
linear 0.180 0.339 0.399 0.403 0.074 0.145
kfobs 0.664 0.728 0.532 0.567 0.163 0.256
kflatent 0.627 0.682 0.537 0.555 0.172 0.267
kflatent_kf 0.648 0.697 0.550 0.572 0.196 0.291
wf 0.466 0.510 0.500 0.536 0.203 0.292
ukf 0.683 0.751 0.565 0.610 0.173 0.277
refh 0.460 0.638 0.667 0.700 0.389 0.436
refh_kf 0.735 0.801 0.702 0.727 0.399 0.459
static rEFH map, training R2: {'x': np.float64(0.654), 'y': np.float64(0.753), 'vx': np.float64(0.776), 'vy': np.float64(0.83), 'ax': np.float64(0.56), 'ay': np.float64(0.636)}
('train', 80), R2 per decoder and variable:
x y vx vy ax ay
linear 0.023 0.064 0.177 0.185 -0.044 0.027
kfobs 0.472 0.462 0.411 0.479 0.074 0.183
kflatent 0.346 0.397 0.395 0.440 0.080 0.156
kflatent_kf 0.527 0.564 0.467 0.542 0.131 0.242
wf 0.248 0.192 0.387 0.387 0.145 0.181
ukf 0.536 0.551 0.426 0.517 0.084 0.217
refh 0.090 0.101 0.427 0.411 0.183 0.217
refh_kf 0.334 0.521 0.518 0.548 0.293 0.320
static rEFH map, training R2: {'x': np.float64(0.327), 'y': np.float64(0.384), 'vx': np.float64(0.47), 'vy': np.float64(0.567), 'ax': np.float64(0.299), 'ay': np.float64(0.347)}
('train', 160), R2 per decoder and variable:
x y vx vy ax ay
linear 0.079 0.194 0.238 0.246 -0.011 0.100
kfobs 0.555 0.629 0.456 0.531 0.097 0.250
kflatent 0.493 0.573 0.450 0.492 0.100 0.219
kflatent_kf 0.483 0.646 0.479 0.541 0.158 0.256
wf 0.415 0.444 0.468 0.489 0.195 0.258
ukf 0.591 0.657 0.479 0.575 0.105 0.277
refh 0.171 0.263 0.532 0.530 0.259 0.301
refh_kf 0.598 0.666 0.615 0.643 0.330 0.393
static rEFH map, training R2: {'x': np.float64(0.423), 'y': np.float64(0.454), 'vx': np.float64(0.647), 'vy': np.float64(0.674), 'ax': np.float64(0.436), 'ay': np.float64(0.474)}
('frac', 0.125), R2 per decoder and variable:
x y vx vy ax ay
linear -0.033 -0.009 0.075 0.047 0.023 0.015
kfobs 0.111 0.151 0.185 0.140 0.065 0.016
kflatent -0.042 -0.007 0.117 0.062 0.014 -0.000
kflatent_kf -0.759 -0.888 -0.039 0.004 -0.037 -0.064
wf 0.072 0.028 0.203 0.185 0.085 0.044
ukf 0.117 0.123 0.188 0.114 0.059 -0.001
refh -0.030 -0.006 0.093 0.044 0.038 0.025
refh_kf 0.123 0.165 0.190 0.145 0.069 0.046
static rEFH map, training R2: {'x': np.float64(0.043), 'y': np.float64(0.043), 'vx': np.float64(0.086), 'vy': np.float64(0.071), 'ax': np.float64(0.069), 'ay': np.float64(0.041)}
('frac', 0.25), R2 per decoder and variable:
x y vx vy ax ay
linear -0.034 0.019 0.087 0.095 0.034 0.041
kfobs -0.012 0.146 0.215 0.299 0.073 0.105
kflatent -0.031 0.074 0.182 0.264 0.047 0.067
kflatent_kf -0.011 0.043 0.181 0.239 0.025 0.045
wf 0.064 0.079 0.265 0.337 0.099 0.132
ukf -0.031 0.169 0.207 0.310 0.064 0.110
refh -0.033 0.006 0.158 0.113 0.079 0.070
refh_kf 0.109 0.128 0.280 0.323 0.143 0.170
static rEFH map, training R2: {'x': np.float64(0.073), 'y': np.float64(0.123), 'vx': np.float64(0.182), 'vy': np.float64(0.183), 'ax': np.float64(0.118), 'ay': np.float64(0.102)}
('frac', 0.5), R2 per decoder and variable:
x y vx vy ax ay
linear 0.010 0.133 0.183 0.197 0.048 0.059
kfobs 0.421 0.587 0.379 0.417 0.116 0.161
kflatent 0.154 0.477 0.332 0.385 0.075 0.105
kflatent_kf 0.033 0.384 0.357 0.341 0.079 0.098
wf 0.241 0.383 0.397 0.407 0.115 0.163
ukf 0.411 0.603 0.380 0.450 0.103 0.175
refh 0.010 0.168 0.390 0.351 0.146 0.142
refh_kf 0.496 0.571 0.492 0.465 0.240 0.266
static rEFH map, training R2: {'x': np.float64(0.158), 'y': np.float64(0.285), 'vx': np.float64(0.488), 'vy': np.float64(0.428), 'ax': np.float64(0.236), 'ay': np.float64(0.231)}
Negative R2 outside 64 ms with 320 s of training¶
fig2-gross-failures is judged at 64 ms with all units and 320 s of training (in
reproduce_makin_2018.ipynb). The other configurations are listed here.
neg = [(k, dec, v, round(r.loc[dec, v], 3)) for k, r in results.items() if k not in (("train", 320), ("frac", 1.0))
for dec in r.index for v in r.columns if r.loc[dec, v] < 0]
print(pd.DataFrame(neg, columns=["config", "decoder", "variable", "R2"]).to_string(index=False))
print("decoders with a negative R2:", sorted({n[1] for n in neg}))
config decoder variable R2
(bin, 16) linear x -0.000
(bin, 16) kflatent_kf ax -0.064
(train, 80) linear ax -0.044
(train, 160) linear ax -0.011
(frac, 0.125) linear x -0.033
(frac, 0.125) linear y -0.009
(frac, 0.125) kflatent x -0.042
(frac, 0.125) kflatent y -0.007
(frac, 0.125) kflatent ay -0.000
(frac, 0.125) kflatent_kf x -0.759
(frac, 0.125) kflatent_kf y -0.888
(frac, 0.125) kflatent_kf vx -0.039
(frac, 0.125) kflatent_kf ax -0.037
(frac, 0.125) kflatent_kf ay -0.064
(frac, 0.125) ukf ay -0.001
(frac, 0.125) refh x -0.030
(frac, 0.125) refh y -0.006
(frac, 0.25) linear x -0.034
(frac, 0.25) kfobs x -0.012
(frac, 0.25) kflatent x -0.031
(frac, 0.25) kflatent_kf x -0.011
(frac, 0.25) ukf x -0.031
(frac, 0.25) refh x -0.033
decoders with a negative R2: ['kflatent', 'kflatent_kf', 'kfobs', 'linear', 'refh', 'ukf']
Outside 64 ms with all units and 320 s of training, negative R2 occurs for linear regression (16 ms, 80 and 160 s, f = 1/8 and 1/4), static KFlatent, KFlatent (+KF) and the static rEFH (f = 1/8 and 1/4), KFobs (f = 1/4, x position), the UKF (f = 1/8, ay, and f = 1/4, x position), and KFlatent (+KF) at 16 ms (x acceleration). Most of these are within 0.05 of zero; the gross failure is KFlatent (+KF) at f = 1/8, whose filter diverges (R2 down to -4.8e7). R2 shown as -0.000 is negative below the printed precision.
fig5-advantage-with-binwidth¶
Reconstruction improves with bin width, and so does the relative advantage of the rEFH (+KF).
OTHERS = ["linear", "kfobs", "kflatent", "kflatent_kf", "wf", "ukf"]
widths = [16, 32, 64, 128]
def advantage(S, v):
best = S.loc[OTHERS, v].max()
return S.loc["refh_kf", v] - best, 100 * (S.loc["refh_kf", v] - best) / best, S.loc[OTHERS, v].idxmax()
rows = []
for w in widths:
S = SNR[("bin", w)]
for v in ("vel", "acc"):
db, pc, who = advantage(S, v)
rows.append(dict(width=w, var=v, refh_kf=S.loc["refh_kf", v], best_other=who, adv_dB=db, adv_pct=pc))
A = pd.DataFrame(rows)
print(A.round(2).to_string(index=False))
def rises(vals, tol):
vals = np.asarray(vals)
return bool(vals[-1] > vals[0] and np.all(np.diff(vals) >= -tol))
verdict = {}
for reading, col, tol in (("a dB", "adv_dB", 0.2), ("b percent", "adv_pct", 2.0)):
ok = True
for v in ("vel", "acc"):
sub = A[A["var"] == v]
snr_ok, adv_ok = rises(sub.refh_kf, 0.2), rises(sub[col], tol)
print(f"reading {reading}, {v}: rEFH (+KF) SNR rises {snr_ok}; advantage rises {adv_ok}")
ok &= snr_ok and adv_ok
verdict[reading] = ok
print("reproduced under:", [k for k, v in verdict.items() if v])
print("SNR of every decoder by width (velocity):"); print(pd.DataFrame({w: SNR[("bin", w)]["vel"] for w in widths}).round(2).to_string())
width var refh_kf best_other adv_dB adv_pct
16 vel 2.96 wf -0.33 -10.14
16 acc 0.95 wf -0.44 -31.67
32 vel 4.02 ukf 0.51 14.63
32 acc 1.69 wf 0.31 22.68
64 vel 5.16 ukf 1.34 35.15
64 acc 2.37 wf 1.04 77.66
128 vel 5.45 ukf 1.59 41.34
128 acc 2.44 wf 1.20 96.08
reading a dB, vel: rEFH (+KF) SNR rises True; advantage rises True
reading a dB, acc: rEFH (+KF) SNR rises True; advantage rises True
reading b percent, vel: rEFH (+KF) SNR rises True; advantage rises True
reading b percent, acc: rEFH (+KF) SNR rises True; advantage rises True
reproduced under: ['a dB', 'b percent']
SNR of every decoder by width (velocity):
16 32 64 128
linear 0.46 0.83 1.39 2.22
kfobs 2.52 3.01 3.30 3.47
kflatent 2.19 2.77 3.19 3.43
kflatent_kf 2.60 3.06 3.40 3.58
wf 3.30 3.28 3.25 3.17
ukf 2.82 3.51 3.82 3.86
refh 0.54 1.68 3.94 5.00
refh_kf 2.96 4.02 5.16 5.45
Result: reproduced. The rEFH (+KF) velocity SNR rises from 3.01 dB at 16 ms to 5.46 dB at 128 ms and acceleration from 0.96 to 2.39 dB, with no fall between neighbouring widths. Its advantage over the best other decoder rises from -0.29 to 1.60 dB for velocity and from -0.44 to 1.14 dB for acceleration (-8.9 to 41.6% and -31.3 to 92.0%), so the claim holds under both readings. At 16 ms the WF beats the rEFH (+KF) for both variables. Not every decoder improves with width: the WF's velocity SNR falls from 3.30 to 3.17 dB.
fig5-derivative-order¶
With increasing derivative order reconstruction quality declines (except that the rEFH is best at velocity), while the relative advantage of the rEFH (+KF) increases.
rows = []
for w in widths:
S = SNR[("bin", w)]
mean_other = S.loc[OTHERS].mean()
pct = {v: advantage(S, v)[1] for v in ("pos", "vel", "acc")}
q = mean_other["pos"] > mean_other["vel"] > mean_other["acc"]
a = pct["pos"] < pct["vel"] < pct["acc"]
rows.append(dict(width=w, mean6_pos=mean_other["pos"], mean6_vel=mean_other["vel"], mean6_acc=mean_other["acc"],
adv_pct_pos=pct["pos"], adv_pct_vel=pct["vel"], adv_pct_acc=pct["acc"], quality_declines=q,
advantage_rises=a, both=q and a, refh_kf_best_var=S.loc["refh_kf"].idxmax()))
D = pd.DataFrame(rows)
print(D.round(2).to_string(index=False))
print("widths where both hold:", int(D.both.sum()), "of 4")
width mean6_pos mean6_vel mean6_acc adv_pct_pos adv_pct_vel adv_pct_acc quality_declines advantage_rises both refh_kf_best_var
16 2.88 2.32 0.51 -5.50 -10.14 -31.67 True False False pos
32 3.39 2.74 0.78 -0.93 14.63 22.68 True True True pos
64 3.79 3.06 0.99 4.80 35.15 77.66 True True True pos
128 4.07 3.29 1.03 15.72 41.34 96.08 True True True pos
widths where both hold: 3 of 4
Result: partially reproduced. The mean SNR of the six other decoders falls from position to velocity to acceleration at all four widths (at 64 ms 3.79, 3.06 and 0.99 dB). The percent advantage of the rEFH (+KF) rises from position to velocity to acceleration at 32, 64 and 128 ms (at 64 ms 4.9, 32.6 and 76.5%), but not at 16 ms, where it is negative and falls (-6.2, -8.9 and -31.3%). Both parts hold at 3 of 4 widths. The claim's exception does not hold here: the rEFH (+KF) decodes position best at every width.
advantage-grows-with-derivative¶
At 64 ms and with each amount of training data, the relative advantage of the rEFH (+KF) over the other decoders increases from position to velocity to acceleration.
for label, res in (("test on all bins after training (main)", results), ("test on bins after 320 s (variant)", results_320)):
rows = []
for tr in (80, 160, 320):
S = summary(M.snr(res[("train", tr)]))
pct = {v: advantage(S, v)[1] for v in ("pos", "vel", "acc")}
rows.append(dict(train_s=tr, adv_pct_pos=pct["pos"], adv_pct_vel=pct["vel"], adv_pct_acc=pct["acc"],
best_other_pos=advantage(S, "pos")[2], rises=pct["pos"] < pct["vel"] < pct["acc"]))
print(label); print(pd.DataFrame(rows).round(2).to_string(index=False))
for tr in (80, 160):
print(f"\nSNR (dB), {tr} s of training, main test set:"); print(summary(M.snr(results[('train', tr)])).round(2).to_string())
test on all bins after training (main)
train_s adv_pct_pos adv_pct_vel adv_pct_acc best_other_pos rises
80 -27.69 8.16 75.34 kflatent_kf True
160 2.13 31.43 74.40 ukf True
320 4.80 35.15 77.66 ukf True
test on bins after 320 s (variant)
train_s adv_pct_pos adv_pct_vel adv_pct_acc best_other_pos rises
80 -29.54 8.03 73.37 ukf True
160 4.41 31.34 74.71 ukf True
320 4.80 35.15 77.66 ukf True
SNR (dB), 80 s of training, main test set:
pos vel acc
linear 0.19 0.87 -0.03
kfobs 2.73 2.56 0.60
kflatent 2.02 2.35 0.55
kflatent_kf 3.43 3.06 0.91
wf 1.08 2.13 0.77
ukf 3.41 2.78 0.72
refh 0.44 2.36 0.97
refh_kf 2.48 3.31 1.59
SNR (dB), 160 s of training, main test set:
pos vel acc
linear 0.65 1.20 0.21
kfobs 3.91 2.97 0.85
kflatent 3.32 2.77 0.77
kflatent_kf 3.69 3.11 1.02
wf 2.44 2.83 1.12
ukf 4.27 3.28 0.95
refh 1.07 3.29 1.43
refh_kf 4.36 4.31 1.95
Result: reproduced. The percent advantage of the rEFH (+KF) rises from position to velocity to acceleration with 80 s (-26.1, 10.9 and 76.4%), 160 s (2.4, 32.5 and 71.3%) and 320 s (4.9, 32.6 and 76.5%) of training. It also rises at each amount when the shorter trainings are scored on the bins after 320 s only.
fig7-affine-position¶
Position SNR of the decoders without a learned kinematic model (linear regression, static KFlatent, WF, static rEFH) is close to an affine function of the fraction of neurons, with correlations 0.999, 0.95, 0.98 and 1.000; static KFlatent is the most curved.
fr = [1 / 8, 1 / 4, 1 / 2, 1.0]
P = pd.DataFrame({f: SNR[("frac", f)]["pos"] for f in fr})
P.columns = ["1/8", "1/4", "1/2", "1"]
print("position SNR (dB), averaged over x and y:"); print(P.round(2).to_string())
static = ["linear", "kflatent", "wf", "refh"]
expected = dict(linear=0.999, kflatent=0.95, wf=0.98, refh=1.0)
out = {}
for reading, xs in (("a f", np.array(fr)), ("b log2 f", np.log2(fr))):
rho = {k: np.corrcoef(xs, P.loc[k].values)[0, 1] for k in static}
within = {k: abs(rho[k] - expected[k]) <= 0.1 * expected[k] for k in static}
lowest = min(rho, key=rho.get)
out[reading] = dict(**{f"rho_{k}": rho[k] for k in static}, all_within=all(within.values()), lowest=lowest,
holds=all(within.values()) and lowest == "kflatent")
print(pd.DataFrame(out).T.to_string())
position SNR (dB), averaged over x and y:
1/8 1/4 1/2 1
linear -0.09 -0.03 0.33 0.80
kfobs 0.61 0.32 3.11 5.00
kflatent -0.10 0.10 1.77 4.08
kflatent_kf -2.61 0.07 1.12 4.26
wf 0.22 0.32 1.65 2.96
ukf 0.56 0.33 3.15 5.63
refh -0.08 -0.06 0.42 1.49
refh_kf 0.68 0.55 3.33 5.90
rho_linear rho_kflatent rho_wf rho_refh all_within lowest holds
a f 0.994636 0.994692 0.986401 0.989682 True wf False
b log2 f 0.955625 0.94797 0.954261 0.910171 True refh False
Result: not reproduced. The correlations of position SNR with f are 0.995 (linear regression), 0.995 (static KFlatent), 0.986 (WF) and 0.991 (static rEFH), and with log2 f 0.956, 0.948, 0.954 and 0.914; all are within 10% of the paper's values. But static KFlatent does not have the lowest correlation under either reading (the WF has with f, the static rEFH with log2 f), so the ordering condition fails. At f = 1/8 and 1/4 (13 and 27 units) the position SNR of these four decoders is within 0.33 dB of zero.
fig7-kinematic-model-saturates¶
Position SNR of the four decoders that learn a kinematic state-transition model (KFobs, KFlatent (+KF), UKF, rEFH (+KF)) rises much faster at first with added neurons, then saturates.
dyn = ["kfobs", "kflatent_kf", "ukf", "refh_kf"]
Pf = pd.DataFrame({f: SNR[("frac", f)]["pos"] for f in fr})
print(Pf.round(2).to_string())
late = (Pf[1.0] - Pf[0.5]) / 0.5
res = {}
for reading, (f0, f1) in (("a 1/8 to 1/4", (1 / 8, 1 / 4)), ("b 1/8 to 1/2", (1 / 8, 1 / 2))):
init = (Pf[f1] - Pf[f0]) / (f1 - f0)
tab = pd.DataFrame({"initial_slope": init, "late_slope": late})
max_static_init = init[static].max()
ok = {k: bool(init[k] > late[k] and init[k] > max_static_init) for k in dyn}
print(f"\nreading {reading} (dB per unit f):"); print(tab.round(2).to_string())
print(f"largest initial slope among the static decoders: {max_static_init:.2f}; per dynamic decoder: {ok}")
res[reading] = all(ok.values())
print("\nreproduced under:", [k for k, v in res.items() if v])
# KFlatent (+KF) without its f = 1/8 value: slope from 1/4 to 1/2, and the f = 1/8 SNR below which
# it passes reading b against its own late slope
k = "kflatent_kf"
print(f"\nKFlatent (+KF): slope 1/4 to 1/2 {(Pf[0.5][k] - Pf[0.25][k]) / 0.25:.2f}, late slope {late[k]:.2f}; "
f"passes reading b against its late slope when its f = 1/8 SNR is below {Pf[0.5][k] - late[k] * 0.375:.2f} dB")
print("f = 1/8 position SNR of the other decoders:", Pf[0.125].drop(k).round(2).to_dict())
0.125 0.250 0.500 1.000
linear -0.09 -0.03 0.33 0.80
kfobs 0.61 0.32 3.11 5.00
kflatent -0.10 0.10 1.77 4.08
kflatent_kf -2.61 0.07 1.12 4.26
wf 0.22 0.32 1.65 2.96
ukf 0.56 0.33 3.15 5.63
refh -0.08 -0.06 0.42 1.49
refh_kf 0.68 0.55 3.33 5.90
reading a 1/8 to 1/4 (dB per unit f):
initial_slope late_slope
linear 0.48 0.93
kfobs -2.35 3.79
kflatent 1.64 4.62
kflatent_kf 21.43 6.27
wf 0.79 2.64
ukf -1.77 4.95
refh 0.15 2.14
refh_kf -1.03 5.15
largest initial slope among the static decoders: 1.64; per dynamic decoder: {'kfobs': False, 'kflatent_kf': True, 'ukf': False, 'refh_kf': False}
reading b 1/8 to 1/2 (dB per unit f):
initial_slope late_slope
linear 1.12 0.93
kfobs 6.67 3.79
kflatent 5.00 4.62
kflatent_kf 9.94 6.27
wf 3.79 2.64
ukf 6.93 4.95
refh 1.32 2.14
refh_kf 7.06 5.15
largest initial slope among the static decoders: 5.00; per dynamic decoder: {'kfobs': True, 'kflatent_kf': True, 'ukf': True, 'refh_kf': True}
reproduced under: ['b 1/8 to 1/2']
KFlatent (+KF): slope 1/4 to 1/2 4.20, late slope 6.27; passes reading b against its late slope when its f = 1/8 SNR is below -1.23 dB
f = 1/8 position SNR of the other decoders: {'linear': -0.09, 'kfobs': 0.61, 'kflatent': -0.1, 'wf': 0.22, 'ukf': 0.56, 'refh': -0.08, 'refh_kf': 0.68}
Result: reproduced, under reading b. From f = 1/8 to 1/2 the position SNR of KFobs, KFlatent (+KF), the UKF and the rEFH (+KF) rises 6.67, 207.33, 6.93 and 6.22 dB per unit f. That is faster than from 1/2 to 1 (3.79, 6.27, 4.95 and 5.79) and faster than any static decoder's initial rise (largest 5.00, static KFlatent). Under reading a (1/8 to 1/4) only KFlatent (+KF) passes: KFobs, the UKF and the rEFH (+KF) decode position slightly worse at 1/4 than at 1/8. KFlatent (+KF) passes under reading b only because its filter diverges at f = 1/8 (position SNR -76.63 dB), a gross failure: from 1/4 to 1/2 its slope is 4.20, below its late slope of 6.27, and it would pass only with an f = 1/8 SNR below -1.23 dB, where every other decoder is between -0.10 and 0.68 dB. Without that value the claim fails under both readings, so the outcome rests on one decoder's failure.
Single-session parts of claims blocked for the other sessions¶
The claims below need a bootstrap over the paper's 49 sessions, which the dandiset does not hold, so they stay blocked. Their point comparisons on this session are printed here; they do not decide an outcome.
fig5-vel-acc-all-bins¶
For velocity and acceleration, the rEFH (+KF) is the best decoder at every bin width (16, 32, 64 and 128 ms).
rows = []
for w in widths:
S = SNR[("bin", w)]
for v in ("vel", "acc"):
rest = S.drop(index="refh_kf")[v]
rows.append(dict(width=w, var=v, refh_kf=S.loc["refh_kf", v], next_best=rest.idxmax(),
margin=S.loc["refh_kf", v] - rest.max()))
print(pd.DataFrame(rows).round(2).to_string(index=False))
width var refh_kf next_best margin
16 vel 2.96 wf -0.33
16 acc 0.95 wf -0.44
32 vel 4.02 ukf 0.51
32 acc 1.69 wf 0.31
64 vel 5.16 refh 1.22
64 acc 2.37 refh 0.51
128 vel 5.45 refh 0.45
128 acc 2.44 refh 0.13
fig5-position-bins¶
For position the leading decoders cannot be told apart at any bin width except 16 ms, where the UKF is better.
rows = []
for w in widths:
S = SNR[("bin", w)]["pos"]
rows.append(dict(width=w, refh_kf=S["refh_kf"], ukf=S["ukf"], kfobs=S["kfobs"], ukf_minus_refh_kf=S["ukf"] - S["refh_kf"],
kfobs_minus_refh_kf=S["kfobs"] - S["refh_kf"], ukf_minus_kfobs=S["ukf"] - S["kfobs"]))
print("position SNR (dB):"); print(pd.DataFrame(rows).round(2).to_string(index=False))
position SNR (dB):
width refh_kf ukf kfobs ukf_minus_refh_kf kfobs_minus_refh_kf ukf_minus_kfobs
16 4.42 4.68 4.13 0.26 -0.29 0.55
32 5.24 5.29 4.70 0.05 -0.55 0.60
64 5.90 5.63 5.00 -0.27 -0.89 0.62
128 6.39 5.52 5.19 -0.87 -1.19 0.33
fig6-training-time¶
The rEFH (+KF) is the best decoder of velocity and acceleration for 80, 160 and 320 s of training data; for position it never differs significantly from the UKF, and beats KFobs at 320 s but loses to it at 80 s.
for label, res in (("test on all bins after training (main)", results), ("test on bins after 320 s (variant)", results_320)):
rows = []
for tr in (80, 160, 320):
S = summary(M.snr(res[("train", tr)]))
row = dict(train_s=tr)
for v in ("vel", "acc"):
rest = S.drop(index="refh_kf")[v]
row[f"{v}_best"] = S[v].idxmax(); row[f"{v}_margin"] = S.loc["refh_kf", v] - rest.max()
row.update(pos_refh_kf=S.loc["refh_kf", "pos"], pos_ukf=S.loc["ukf", "pos"], pos_kfobs=S.loc["kfobs", "pos"],
pos_minus_ukf=S.loc["refh_kf", "pos"] - S.loc["ukf", "pos"], pos_minus_kfobs=S.loc["refh_kf", "pos"] - S.loc["kfobs", "pos"])
rows.append(row)
print(label); print(pd.DataFrame(rows).round(2).to_string(index=False))
test on all bins after training (main)
train_s vel_best vel_margin acc_best acc_margin pos_refh_kf pos_ukf pos_kfobs pos_minus_ukf pos_minus_kfobs
80 refh_kf 0.25 refh_kf 0.62 2.48 3.41 2.73 -0.93 -0.25
160 refh_kf 1.02 refh_kf 0.52 4.36 4.27 3.91 0.09 0.45
320 refh_kf 1.22 refh_kf 0.51 5.90 5.63 5.00 0.27 0.89
test on bins after 320 s (variant)
train_s vel_best vel_margin acc_best acc_margin pos_refh_kf pos_ukf pos_kfobs pos_minus_ukf pos_minus_kfobs
80 refh_kf 0.25 refh_kf 0.63 2.41 3.43 2.70 -1.01 -0.28
160 refh_kf 0.98 refh_kf 0.51 4.60 4.41 4.02 0.19 0.58
320 refh_kf 1.22 refh_kf 0.51 5.90 5.63 5.00 0.27 0.89
fig7-vel-acc-neurons¶
For velocity and acceleration, using a fraction 1/4, 1/2 or all of the neurons, the rEFH (+KF) outperforms all its rivals, and its percent gain in SNR grows with the number of neurons and with derivative order.
Percent gain under reading (a), over the best decoder other than the two rEFH decoders, and reading (b), over each such decoder.
rows, pct_each = [], {}
for f in (1 / 4, 1 / 2, 1.0):
S = SNR[("frac", f)]
for v in ("vel", "acc"):
rest = S.drop(index="refh_kf")[v]
best = S.loc[OTHERS, v].max()
rows.append(dict(f=f, var=v, refh_kf=S.loc["refh_kf", v], next_best=rest.idxmax(), margin=S.loc["refh_kf", v] - rest.max(),
pct_over_best_other=100 * (S.loc["refh_kf", v] - best) / best))
pct_each[(f, v)] = 100 * (S.loc["refh_kf", v] - S.loc[OTHERS, v]) / S.loc[OTHERS, v]
print(pd.DataFrame(rows).round(2).to_string(index=False))
print("\npercent gain over each decoder (reading b):"); print(pd.DataFrame(pct_each).round(1).to_string())
f var refh_kf next_best margin pct_over_best_other
0.25 vel 1.56 wf 0.00 0.16
0.25 acc 0.74 wf 0.20 38.24
0.50 vel 2.83 ukf 0.49 21.07
0.50 acc 1.27 refh 0.59 93.71
1.00 vel 5.16 refh 1.22 35.15
1.00 acc 2.37 refh 0.51 77.66
percent gain over each decoder (reading b):
0.25 0.50 1.00
vel acc vel acc vel acc
linear 278.1 342.4 209.2 428.8 269.7 508.2
kfobs 20.5 82.5 28.2 95.5 56.1 133.4
kflatent 41.9 189.3 46.4 207.4 61.6 141.9
kflatent_kf 52.1 377.5 51.6 215.0 51.7 120.7
wf 0.2 38.2 26.5 94.2 58.7 77.7
ukf 19.4 85.6 21.1 93.7 35.2 113.7
fig7-position-neurons¶
For position the rEFH's advantage is significant at p < 0.05 at f = 1/4 and 1/2, but not at f = 1, where the UKF comes close, nor at f = 1/8, where all decoders collapse.
rows = []
for f in fr:
S = SNR[("frac", f)]["pos"]
rest = S.drop("refh_kf")
rows.append(dict(f=f, refh_kf=S["refh_kf"], next_best=rest.idxmax(), margin=S["refh_kf"] - rest.max(),
minus_ukf=S["refh_kf"] - S["ukf"]))
print("position SNR (dB):"); print(pd.DataFrame(rows).round(2).to_string(index=False))
position SNR (dB): f refh_kf next_best margin minus_ukf 0.12 0.68 kfobs 0.07 0.12 0.25 0.55 ukf 0.21 0.21 0.50 3.33 ukf 0.17 0.17 1.00 5.90 ukf 0.27 0.27
Result: blocked, single-session values only. These claims need the bootstrap over the paper's 49 sessions. On this session:
- fig5-vel-acc-all-bins: the rEFH (+KF) has the highest velocity and acceleration SNR at 32, 64 and 128 ms (margins 0.53 and 0.32, 1.19 and 0.45, 0.45 and 0.12 dB), but at 16 ms the WF is higher for both (3.30 against 3.01 dB, and 1.39 against 0.96 dB).
- fig5-position-bins: at 16 ms the UKF's position SNR exceeds the rEFH (+KF)'s by 0.29 dB, the direction the paper reports. At 32, 64 and 128 ms the differences between the rEFH (+KF), the UKF and KFobs are 0.01 to 1.22 dB in size; whether any is significant needs the bootstrap.
- fig6-training-time: part (i) holds on point values, the rEFH (+KF) having the highest velocity and acceleration SNR at 80, 160 and 320 s (velocity margins 0.33, 1.01 and 1.19 dB, acceleration 0.68, 0.47 and 0.45 dB). For position it is 0.87 dB below the UKF at 80 s and 0.10 and 0.28 dB above it at 160 and 320 s, and against KFobs it is above at 320 s (0.90 dB) and below at 80 s (0.20 dB), the directions part (iii) states. Scoring on the bins after 320 s gives the same signs.
- fig7-vel-acc-neurons: at f = 1/4 the rEFH (+KF) ties the WF for velocity (1.56 dB each, 0.00 dB apart) and exceeds it by 0.20 dB for acceleration; at f = 1/2 and 1 it is the best decoder of both. Its percent gain over the best decoder other than the two rEFH decoders (reading a) rises with f for velocity (0.2, 19.1 and 32.6%) but not for acceleration (38.2, 87.1 and 76.5%); it is larger for acceleration than velocity at each f. Over each decoder (reading b) the gain also falls somewhere, for example over the UKF for velocity from 19.4% at f = 1/4 to 19.1% at f = 1/2.
- fig7-position-neurons: the rEFH (+KF) has the highest position SNR at f = 1/8 (by 0.07 dB; every decoder but the diverged KFlatent (+KF) within 0.7 dB of zero), at f = 1/4 (by 0.21 dB) and at f = 1 (0.28 dB over the UKF), and is 0.15 dB below the UKF at f = 1/2, against the paper's advantage there.
Summary¶
| Claim | Outcome |
|---|---|
| fig5-advantage-with-binwidth | reproduced |
| fig5-derivative-order | partially reproduced (3 of 4 widths; not at 16 ms) |
| advantage-grows-with-derivative | reproduced |
| fig7-affine-position | not reproduced (correlations within tolerance; static KFlatent not the most curved) |
| fig7-kinematic-model-saturates | reproduced under reading b, only because KFlatent (+KF) diverges at f = 1/8 |
| fig5-vel-acc-all-bins, fig5-position-bins, fig6-training-time, fig7-vel-acc-neurons, fig7-position-neurons | blocked: bootstrap across the paper's 49 sessions, which the dandiset does not hold (single-session values above) |
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 (dated 2017-02-02, not in the paper's
session list) in place of 49, so each mean over sessions is this session's value. Sorted units
only, 111 above 0.5 Hz, without hash units; a fraction f keeps the first floor(111 f) of them
(13, 27 and 55 units). For 80 and 160 s of training the test bins are all bins after the
training bins, as in the paper's code; the paper does not say. The decoders are the same
re-implementations as in reproduce_makin_2018.ipynb; KFlatent's EM reached the 100-epoch limit
in every configuration without meeting the 1/100 of a percent criterion (printed with the fits).
What was found about the dataset. As in reproduce_makin_2018.ipynb: one session, no hash
units, three undocumented 0.1 s gaps, and a test file without kinematics. The same 111 units pass
the rate filter at every bin width.
print(f"notebook run time {(time.time() - T0) / 60:.1f} min")
notebook run time 10.1 min