v0.3.0: Full Matlab parity — methods, examples, and decoding algorithms

README.md

@@ -12,8 +12,9 @@ One of nSTAT's key strengths is point-process generalized linear models for spik
 train signals that provide a formal statistical framework for processing signals
 recorded from ensembles of single neurons. It also has extensive support for model
 fitting, model-order analysis, and adaptive decoding — including state-space GLM
-(SSGLM) estimation via EM, unscented Kalman filtering (UKF), and hybrid
-discrete/continuous point-process filters.
+(SSGLM) estimation via EM, unscented Kalman filtering (UKF), goal-directed
+point-process adaptive filters (PPAF), and hybrid discrete/continuous
+point-process filters (PPHF).
 
 Although created with neural signal processing in mind, nSTAT can be used as a
 generic tool for analyzing any types of discrete and continuous signals, and thus

docs/ClassDefinitions.md

@@ -102,12 +102,21 @@ Primary notebook: `notebooks/nstCollExamples.ipynb`
 **Time operations**:
 `shiftTime`, `setMinTime`, `setMaxTime`
 
+**Data export**:
+`dataToMatrix`, `resample`
+
+**PSTH**:
+`psth`, `psthGLM`, `estimateVarianceAcrossTrials`, `psthBars`
+
 **SSGLM (state-space GLM)**:
 `ssglm`, `ssglmFB`
 
 **Basis generation**:
 `generateUnitImpulseBasis`
 
+**Plotting**:
+`plot`
+
 ### `History` (`nstat.History`)
 
 Primary notebook: `notebooks/HistoryExamples.ipynb`
@@ -149,7 +158,8 @@ Primary notebook: `notebooks/TrialExamples.ipynb`
 `flattenMask`
 
 **Utilities**:
-`shiftCovariates`, `resampleEnsColl`, `restoreToOriginal`, `plot`
+`shiftCovariates`, `resampleEnsColl`, `restoreToOriginal`, `getAllLabels`,
+`plot`
 
 ### `TrialConfig` (`nstat.TrialConfig`)
 
@@ -230,19 +240,23 @@ Alias of `FitSummary`. Aggregates multiple `FitResult` objects.
 
 **Plotting**:
 `plotIC`, `plotAIC`, `plotBIC`, `plotlogLL`, `plotResidualSummary`,
-`plotSummary`, `boxPlot`
+`plotSummary`, `boxPlot`, `plotAllCoeffs`, `plot3dCoeffSummary`,
+`plot2dCoeffSummary`, `plotKSSummary`
 
 ### `DecodingAlgorithms` (`nstat.DecodingAlgorithms`)
 
 Primary notebook: `notebooks/DecodingExample.ipynb`
 
 **Point-process decode filters**:
 `PPDecodeFilterLinear`, `PPDecodeFilter`, `PPHybridFilterLinear`,
-`ComputeStimulusCIs`
+`ComputeStimulusCIs`, `computeSpikeRateCIs`
 
 **Kalman and unscented Kalman filters**:
 `kalman_filter`, `PP_fixedIntervalSmoother`, `ukf`
 
+**Helper methods**:
+`PPDecode_predict`, `PPDecode_updateLinear`
+
 **State-space GLM (SSGLM) — EM algorithm**:
 `PPSS_EStep`, `PPSS_MStep`, `PPSS_EM`, `PPSS_EMFB`
 

docs/PaperOverview.md

@@ -88,10 +88,16 @@ map to:
 **Point-process adaptive filters (Section 2.5)**:
 
 - `DecodingAlgorithms.PPDecodeFilterLinear` — linear-CIF point-process
-  adaptive filter for continuous stimulus decoding.
+  adaptive filter for continuous stimulus decoding. Supports goal-directed
+  decoding via backward information filter (Srinivasan et al. 2006) when
+  `yT` and `PiT` target parameters are provided.
 - `DecodingAlgorithms.PPDecodeFilter` — general CIF version using symbolic
   gradients/Jacobians.
-- `DecodingAlgorithms.ComputeStimulusCIs` — stimulus confidence intervals.
+- `DecodingAlgorithms.ComputeStimulusCIs` — stimulus confidence intervals
+  via Monte Carlo sampling (dual-path: 4-D SSGLM cross-trial + 3-D smoother
+  z-score).
+- `DecodingAlgorithms.computeSpikeRateCIs` — spike rate confidence intervals
+  and pairwise significance testing across trials.
 - `DecodingAlgorithms.PP_fixedIntervalSmoother` — fixed-interval smoother
   for off-line smoothing of decode estimates.
 
@@ -100,6 +106,8 @@ map to:
 - `DecodingAlgorithms.PPHybridFilterLinear` — joint discrete/continuous state
   estimation combining point-process observations with a hidden Markov model
   over discrete states and Kalman filtering over continuous kinematics.
+  Supports goal-directed decoding via per-model backward information filters
+  when `yT` and `PiT` are provided.
 
 **Kalman and UKF filters**:
 

examples/README.md

@@ -1,17 +1,44 @@
-# Python nSTAT Examples
+# nSTAT Python Examples
 
-## Basic examples
+## Paper Examples (Self-Contained)
+
+Five self-contained scripts mirroring the MATLAB paper examples. Each
+generates publication-quality figures and supports `--export-figures`.
 
 ```bash
-python3 examples/basic_data_workflow.py
-python3 examples/fit_poisson_glm.py
-python3 examples/simulate_population_psth.py
+python examples/paper/example01_mepsc_poisson.py --export-figures
+python examples/paper/example02_whisker_stimulus_thalamus.py --export-figures
+python examples/paper/example03_psth_and_ssglm.py --export-figures
+python examples/paper/example04_place_cells_continuous_stimulus.py --export-figures
+python examples/paper/example05_decoding_ppaf_pphf.py --export-figures
 ```
 
-## Paper-style example workflow
+| Example | Focus | Paper Section |
+|---|---|---|
+| 01 | mEPSC Poisson models (constant vs piecewise baseline) | 2.3.1 |
+| 02 | Whisker stimulus GLM with lag and history selection | 2.3.2 |
+| 03 | PSTH and SSGLM across-trial dynamics | 2.3.3-2.3.4 |
+| 04 | Place-cell receptive fields (Gaussian vs Zernike) | 2.3.5 |
+| 05 | PPAF and hybrid filter decoding | 2.5-2.6 |
+
+## Basic Examples
 
 ```bash
-python3 examples/nstat_paper_examples.py --repo-root ..
+python examples/basic_data_workflow.py
+python examples/fit_poisson_glm.py
+python examples/simulate_population_psth.py
 ```
 
-This mirrors key analyses described in the nSTAT paper using the bundled Python APIs.
+## README Examples (Quick Checks)
+
+```bash
+python examples/readme_examples/example1_multitaper_and_spectrogram.py
+python examples/readme_examples/example2_simulate_cif_spiketrain_10s.py
+python examples/readme_examples/example3_nstcoll_raster_from_example2.py
+```
+
+## Jupyter Notebooks
+
+All 29 class-tutorial and data-analysis notebooks are in `notebooks/`.
+They mirror the MATLAB helpfile examples one-to-one. See
+[docs/Examples.md](../docs/Examples.md) for the full index.

examples/paper/example01_mepsc_poisson.py

@@ -128,33 +128,14 @@ def run_example01(*, export_figures: bool = False, export_dir: Path | None = Non
     print(f"  AIC: {resultConst.AIC}")
     print(f"  BIC: {resultConst.BIC}")
 
-    # --- Figure 1: Constant Mg2+ raster + diagnostics + lambda ---
-    fig1, axes1 = plt.subplots(2, 2, figsize=(14, 9))
-
-    # Subplot 1: Neural raster
-    ax = axes1[0, 0]
-    spikeCollConst.plot(handle=ax)
-    ax.set_title("Neural Raster with constant Mg$^{2+}$ Concentration",
-                 fontweight="bold", fontsize=12)
-    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
-    ax.set_ylabel("mEPSCs", fontsize=12, fontweight="bold")
-    ax.set_yticks([0, 1])
-
-    # Subplot 2: Inverse Gaussian Transform (autocorrelation of rescaled residuals)
-    ax = axes1[0, 1]
-    resultConst.plotInvGausTrans(handle=ax)
-
-    # Subplot 3: KS plot
-    ax = axes1[1, 0]
-    resultConst.KSPlot(handle=ax)
-
-    # Subplot 4: Lambda estimate
-    ax = axes1[1, 1]
-    resultConst.lambda_signal.plot(handle=ax)
-    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
-    ax.legend(["$\\lambda_{const}$"], loc="upper right")
-
-    fig1.suptitle("Example 01 — Figure 1: Constant Mg$^{2+}$ Summary", fontsize=14, fontweight="bold")
+    # --- Figure 1: Constant Mg2+ diagnostics (Matlab-matching plotResults) ---
+    # Matlab calls resultConst.plotResults which creates a 2x4 grid:
+    #   [KSPlot (2-wide)]  [InvGausTrans]  [SeqCorr]
+    #   [plotCoeffs (2-wide)]               [plotResidual (2-wide)]
+    fig1 = plt.figure(figsize=(14, 9))
+    resultConst.plotResults(handle=fig1)
+    fig1.suptitle("Example 01 — Figure 1: Constant Mg$^{2+}$ Summary",
+                  fontsize=14, fontweight="bold")
     fig1.tight_layout()
     figure_files.extend(_maybe_export(fig1, export_dir, "fig01_constant_mg_summary"))
 
@@ -198,8 +179,12 @@ def run_example01(*, export_figures: bool = False, export_dir: Path | None = Non
     print("\n=== Part 3: Piecewise Baseline Model Comparison ===")
 
     # Build piecewise indicator covariates
-    timeInd1 = np.searchsorted(timeWashout, 495.0)
-    timeInd2 = np.searchsorted(timeWashout, 765.0)
+    # Matlab: find(time<495,1,'last') — last index strictly before 495
+    # np.searchsorted gives first index >= 495, so subtract 1 isn't needed
+    # because Python slice [:idx] is exclusive. But Matlab's 1:timeInd1 is
+    # inclusive, so we need searchsorted(..., side='right') to include 495.
+    timeInd1 = np.searchsorted(timeWashout, 495.0, side="right")
+    timeInd2 = np.searchsorted(timeWashout, 765.0, side="right")
     N = len(timeWashout)
 
     constantRate = np.ones((N, 1))
@@ -233,34 +218,12 @@ def run_example01(*, export_figures: bool = False, export_dir: Path | None = Non
     print(f"  AIC: {resultWashout.AIC}")
     print(f"  BIC: {resultWashout.BIC}")
 
-    # --- Figure 3: Piecewise model diagnostics + lambda comparison ---
-    fig3, axes3 = plt.subplots(2, 2, figsize=(14, 9))
-
-    # Subplot 1: Raster with epoch boundaries
-    ax = axes3[0, 0]
-    spikeCollWashout.plot(handle=ax)
-    ax.set_title("Neural Raster with decreasing Mg$^{2+}$ Concentration",
-                 fontweight="bold", fontsize=12)
-    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
-    ax.axvline(495.0, color="r", linewidth=4)
-    ax.axvline(765.0, color="r", linewidth=4)
-
-    # Subplot 2: Inverse Gaussian Transform
-    ax = axes3[0, 1]
-    resultWashout.plotInvGausTrans(handle=ax)
-
-    # Subplot 3: KS plot
-    ax = axes3[1, 0]
-    resultWashout.KSPlot(handle=ax)
-
-    # Subplot 4: Lambda comparison
-    ax = axes3[1, 1]
-    resultWashout.lambda_signal.plot(handle=ax)
-    ax.set_ylim(0, 5)
-    ax.set_xlabel("time [s]", fontsize=12, fontweight="bold")
-    ax.legend(["$\\lambda_{const}$", "$\\lambda_{const-epoch}$"], loc="upper right")
-
-    fig3.suptitle("Example 01 — Figure 3: Piecewise Baseline Comparison", fontsize=14, fontweight="bold")
+    # --- Figure 3: Piecewise model diagnostics (Matlab-matching plotResults) ---
+    # Matlab calls resultWashout.plotResults which creates the same 2x4 grid.
+    fig3 = plt.figure(figsize=(14, 9))
+    resultWashout.plotResults(handle=fig3)
+    fig3.suptitle("Example 01 — Figure 3: Piecewise Baseline Comparison",
+                  fontsize=14, fontweight="bold")
     fig3.tight_layout()
     figure_files.extend(_maybe_export(fig3, export_dir, "fig03_piecewise_baseline_comparison"))