nstat.DecodingAlgorithms

class nstat.DecodingAlgorithms[source]

Bases: object

Static-method library for neural decoding and state-space estimation.

Provides Kalman filtering/smoothing, point-process adaptive filtering (PPAF), hybrid discrete–continuous decoding, unscented Kalman filtering (UKF), state-space GLM EM algorithms (SSGLM), and mixed point-process / continuous-observation EM algorithms (PPLFP_*, with the deprecated mPPCO_* aliases).

All methods are @staticmethod; no instance is required. Method signatures follow the Matlab DecodingAlgorithms class as closely as possible.

See also

CIF

Conditional intensity function objects used by PPDecodeFilter.

Analysis

High-level fitting routines that call these decoders.

__init__()
static PPDecodeFilter(A, Q, Px0, dN, lambdaCIFColl, binwidth=0.001, x0=None, Pi0=None, yT=None, PiT=None, estimateTarget=0, Wconv=None)[source]

Point-process adaptive filter using CIF object collection.

Runs the full forward filter loop, evaluating CIF objects at each time step. When yT / PiT are supplied, delegates to the linear-parameter variant with goal-state estimation.

Returns (x_p, W_p, x_u, W_u, xT, WT, MT, MT_cov).

static PPDecodeFilterLinear(*args, **kwargs)[source]

Point-process adaptive filter using linear GLM parameters.

Dispatches to _ppdecode_filter_linear when a fitType string is present, otherwise falls back to kalman_filter. Matches the Matlab DecodingAlgorithms.PPDecodeFilterLinear signature.

static PP_fixedIntervalSmoother(A, Q, dN, lags, mu, beta, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, x0=None, Pi0=None)[source]

Point-process fixed-interval (fixed-lag) smoother.

Runs PPDecode_updateLinear forward, then applies backward RTS-style smoothing at each step with the specified number of lags. Returns (x_pLag, W_pLag, x_uLag, W_uLag).

static PPDecode_predict(x_u, W_u, A, Q, Wconv=None)[source]

Point-process adaptive filter predict step.

Returns (x_p, W_p) — predicted state and covariance.

static PPDecode_update(x_p, W_p, dN, lambdaIn, binwidth=0.001, time_index=0, WuConv=None)[source]

Point-process adaptive filter update step using CIF objects.

Evaluates symbolic CIF gradients and Jacobians for the Newton-step posterior update. Returns (x_u, W_u, lambda_delta).

static PPDecode_updateLinear(x_p, W_p, dN, mu, beta, fitType='poisson', gamma=None, HkAll=None, time_index=0, WuConv=None)[source]

Point-process adaptive filter update step using linear parameters.

Uses mu, beta, and optional gamma history coefficients instead of CIF objects. Returns (x_u, W_u, lambda_delta). As MATLAB’s PPDecode_updateLinear, lambda_delta is exp(linTerm) (binomial exp./(1+exp)), unclipped, with NaN / Inf entries set to 1 (this port clipped linTerm to [-20, 20]).

static PPHybridFilterLinear(A, Q, p_ij, Mu0, dN, mu, beta, fitType='poisson', binwidth=0.001, gamma=None, windowTimes=None, x0=None, Pi0=None, yT=None, PiT=None, estimateTarget=0, MinClassificationError=0)[source]

Hybrid point-process filter with discrete-mode switching (linear parameters).

Combines multiple linear state-space models (one per discrete mode) with a Markov transition matrix p_ij. At each time step, per-mode PPAF updates are merged using posterior mode probabilities.

Returns (x_p, W_p, x_u, W_u, xT, WT, S_est, X_est, W_est).

static PPHybridFilter(A, Q, p_ij, Mu0, dN, lambdaCIFColl, binwidth=0.001, x0=None, Pi0=None, yT=None, PiT=None, estimateTarget=0, MinClassificationError=0)[source]

Hybrid point-process filter with CIF-object evaluation.

Unlike PPHybridFilterLinear() which takes pre-extracted linear parameters (mu, beta, gamma), this method evaluates CIF objects directly via their evalLambdaDelta / evalGradient* / evalJacobian* methods. This supports nonlinear conditional intensity specifications.

Falls back to the linear path when the target-estimation branch is active (yT / PiT / estimateTarget supplied), matching Matlab behaviour.

static ukf(fstate, x, P, hmeas, z, Q, R)[source]

Unscented Kalman Filter for nonlinear systems.

One-step update matching Matlab DecodingAlgorithms.ukf.

System model (additive noise):

x_{k+1} = fstate(x_k) + w_k,   w ~ N(0, Q)
z_k     = hmeas(x_k)  + v_k,   v ~ N(0, R)
Parameters:

  • fstate (callable (L,) -> (L,))

  • x ((L,) prior state estimate)

  • P ((L, L) prior covariance)

  • hmeas (callable (L,) -> (m,))

  • z ((m,) measurement)

  • Q ((L, L) process noise covariance)

  • R ((m, m) measurement noise covariance)

Returns:

  • x ((L,) posterior state estimate)

  • P ((L, L) posterior covariance)

static ukf_ut(f, X, Wm, Wc, n, R)[source]

Unscented transformation.

Parameters:

  • f (callable mapping (L,) -> (n,))

  • X ((L, 2L+1) sigma points)

  • Wm ((2L+1,) weights)

  • Wc ((2L+1,) weights)

  • n (output dimensionality)

  • R ((n, n) additive covariance)

Returns:

  • y ((n,) transformed mean)

  • Y ((n, 2L+1) transformed sigma points)

  • P ((n, n) transformed covariance)

  • Y1 ((n, 2L+1) deviations)

static ukf_sigmas(x, P, c)[source]

Generate sigma points around reference point x.

Parameters:

  • x ((L,) state vector)

  • P ((L, L) covariance)

  • c (scaling coefficient)

Returns:

X

Return type:

(L, 2L+1) sigma-point matrix

static linear_decode(spike_counts, stimulus)[source]

Ordinary-least-squares linear decoder (spike counts → stimulus).

Returns a dict with keys 'coefficients', 'decoded', 'residual', and 'ci' (95 % confidence band).

Parameters:
Return type:

dict[str, ndarray]

static kalman_filter(A=None, C=None, Pv=None, Pw=None, Px0=None, x0=None, y=None, GnConv=None, *, observations=None, transition=None, observation_matrix=None, q_cov=None, r_cov=None, p0=None)[source]

Discrete-time Kalman filter.

Accepts both the MATLAB-compatible signature:

kalman_filter(A, C, Pv, Pw, Px0, x0, y, GnConv)

and the Pythonic keyword signature:

kalman_filter(observations=y, transition=A,
              observation_matrix=C, q_cov=Pv,
              r_cov=Pw, x0=x0, p0=Px0)

When MATLAB-style positional args are provided, delegates to _kalman_filter_matlab. Otherwise runs a simple Pythonic Kalman filter returning dict(state=..., cov=...).

Returns:

state : (N, Dx) — updated (posterior) state estimates. cov : (N, Dx, Dx) — updated covariances.

Return type:

dict with keys

static kalman_predict(x_u, Pe_u, A, Pv, GnConv=None)[source]

Kalman filter predict step: x_p = A x_u, Pe_p = A Pe A' + Pv.

When the converged-gain GnConv is provided the predict-step covariance is taken from the previously-filtered Pe_u (which is the steady-state covariance under the converged gain) rather than recomputed. Audit finding H1 (Agent 5): the prior implementation used GnConv itself as the covariance — but GnConv is a gain matrix (shape Dx x Dy), not a covariance (Dx x Dx). The _kalman_filter_matlab path at lines 611-614 had this right; the standalone kalman_predict diverged. Matches MATLAB KalmanFilter.m:91 (Pe_p = Pe_u when Gn_conv set).

static kalman_update(x_p, Pe_p, C, Pw, y, GnConv=None)[source]

Kalman filter update step: incorporate observation y and return (x_u, Pe_u, G).

static kalman_smootherFromFiltered(A, x_p, Pe_p, x_u, Pe_u)[source]

RTS backward smoother from precomputed filter estimates.

Returns (x_N, P_N, Ln) — smoothed states, covariances, and backward Kalman gains.

static kalman_smoother(A, C, Pv, Pw, Px0, x0, y)[source]

Run a Kalman filter followed by an RTS smoother.

Returns (x_N, P_N, Ln, x_p, Pe_p, x_u, Pe_u).

static kalman_fixedIntervalSmoother(A, C, Pv, Pw, Px0, x0, y, lags)[source]

Kalman fixed-interval (fixed-lag) smoother via augmented state.

Matches the Matlab implementation: builds an augmented state of dimension (1 + lags) * n_x and runs kalman_filter on the augmented system. The lagged portion of the augmented state gives the exact smoothed estimate at lag lags.

Parameters:

  • A (arrays) – System matrices (may be time-varying 3-D arrays).

  • C (arrays) – System matrices (may be time-varying 3-D arrays).

  • Pv (arrays) – System matrices (may be time-varying 3-D arrays).

  • Pw (arrays) – System matrices (may be time-varying 3-D arrays).

  • Px0 ((Dx, Dx) — initial covariance.)

  • x0 ((Dx,) — initial state.)

  • y ((N, Dy) or (Dy, N) — observations (auto-detected layout).)

  • lags (int — number of smoothing lags.)

Returns:

  • x_pLag ((N+1, Dx) — predicted states at the lagged component (time-major).)

  • Pe_pLag ((N+1, Dx, Dx) — predicted covariances at the lagged component.)

  • x_uLag ((N, Dx) — updated states at the lagged component.)

  • Pe_uLag ((N, Dx, Dx) — updated covariances at the lagged component.)

static PPSS_EMFB(A, Q0, x0, dN, fitType, delta, gamma0, windowTimes, numBasis, neuronName=None)[source]

EM Forward-Backward algorithm for state-space GLM.

Wraps PPSS_EM in a forward-backward-forward cycle for improved convergence.

Parameters:

  • A ((R, R) state transition matrix (typically identity for random walk))

  • Q0 ((R,) initial state noise variance)

  • x0 ((R,) initial state coefficients)

  • dN ((K, N) binary spike observations (K trials, N time bins))

  • fitType ('poisson' or 'binomial')

  • delta (time bin width in seconds)

  • gamma0 ((J,) initial history coefficients)

  • windowTimes (array of history window boundaries)

  • numBasis (number of basis functions R)

  • neuronName (identifier for the neuron (for labeling))

Returns:

  • xKFinal ((R, K) estimated state trajectories)

  • WKFinal ((R, R, K) estimated state covariances)

  • WkuFinal ((R, R, K, K) cross-trial covariances)

  • Qhat ((R,) estimated state noise variance)

  • gammahat ((J,) estimated history coefficients)

  • fitResults (FitResult object with goodness-of-fit diagnostics)

  • stimulus ((R, K) estimated stimulus effect)

  • stimCIs ((R, K, 2) stimulus confidence intervals)

  • logll (float, log-likelihood at convergence)

  • QhatAll (parameter history)

  • gammahatAll (parameter history)

  • nIter (total EM iterations)

static PPSS_EM(A, Q0, x0, dN, fitType, delta, gamma0, windowTimes, numBasis, HkAll)[source]

Inner EM loop for state-space GLM.

Parameters:

  • A ((R, R) state transition)

  • Q0 ((R,) initial state noise variance)

  • x0 ((R,) initial state)

  • dN ((K, N) observations)

  • fitType ('poisson' or 'binomial')

  • delta (time bin width)

  • gamma0 ((J,) initial history coefficients)

  • windowTimes (history window boundaries)

  • numBasis (number of basis functions)

  • HkAll (precomputed history matrices)

Return type:

xKFinal, WKFinal, WkuFinal, Qhat, gammahat, logll, QhatAll, gammahatAll, nIter, negLL

static PPSS_EStep(A, Q, x0, dN, HkAll, fitType, delta, gamma, numBasis)[source]

E-step: Forward Kalman filter + backward RTS smoother + cross-covariance.

Parameters:

  • A ((R, R) state transition matrix)

  • Q ((R,) or (R, R) state noise covariance (diagonal vector or matrix))

  • x0 ((R,) initial state)

  • dN ((K, N) binary spike observations (K trials, N time bins))

  • HkAll (list of K arrays, each (N, J) history design matrices)

  • fitType ('poisson' or 'binomial')

  • delta (time bin width)

  • gamma ((J,) history coefficients)

  • numBasis (number of basis functions R)

Return type:

x_K, W_K, Wku, logll, sumXkTerms, sumPPll

static PPSS_MStep(dN, HkAll, fitType, x_K, W_K, gamma, delta, sumXkTerms, windowTimes)[source]

M-step: Update Q via closed form, gamma via Newton-Raphson.

Parameters:

  • dN ((K, N))

  • HkAll (list of K arrays (N, J))

  • fitType ('poisson' or 'binomial')

  • x_K ((R, K) smoothed states)

  • W_K ((R, R, K) smoothed covariances)

  • gamma ((J,) current history coefficients)

  • delta (time bin width)

  • sumXkTerms ((R, R) sufficient statistics from E-step)

  • windowTimes (array of history window boundaries)

Returns:

  • Qhat ((R,) updated state noise variance (diagonal))

  • gamma_new ((J,) updated history coefficients)

static prepareEMResults(fitType, neuronNumber, dN, HkAll, xK, WK, Q, gamma, windowTimes, delta, informationMatrix, logll)[source]

Package SSGLM EM results into a FitResult object.

static ComputeStimulusCIs(fitType, xK, Wku, delta, Mc=None, alphaVal=0.05)[source]

Confidence intervals for stimulus estimate.

When Wku is a 4-D array (numBasis, numBasis, K, K) (SSGLM cross-trial covariance), uses Monte Carlo sampling matching Matlab’s DecodingAlgorithms.ComputeStimulusCIs.

When Wku is a 3-D array (N, Dx, Dx) (e.g. smoother output), falls back to z-score Gaussian CIs with inverse-link transform.

Parameters:

  • fitType (str) – 'poisson' or 'binomial'.

  • xK (array) – Smoothed state estimates. Shape (numBasis, K) for SSGLM or (N, Dx) for smoother output.

  • Wku (array) – Covariance. Shape (numBasis, numBasis, K, K) for SSGLM or (N, Dx, Dx) for smoother output.

  • delta (float) – Time-step size in seconds.

  • Mc (int, optional) – Number of Monte Carlo draws (default 3000). Ignored when using z-score fallback.

  • alphaVal (float, optional) – Significance level for two-sided CIs (default 0.05).

Returns:

  • CIs (array, shape (..., 2)) – Lower and upper confidence bounds.

  • stimulus (array) – Point estimate (inverse-link of xK, divided by delta for MC mode).

static estimateInfoMat(fitType, dN, HkAll, A, x0, xK, WK, Wku, Q, gamma, windowTimes, SumXkTerms, delta, Mc=500)[source]

Observed information matrix via Louis’ identity with Monte Carlo.

Computes I_obs = I_complete - I_missing where I_missing is estimated by MC sampling from the smoothing distribution.

static computeSpikeRateCIs(xK, Wku, dN, t0, tf, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, Mc=500, alphaVal=0.05)[source]

Monte Carlo spike-rate confidence intervals.

Computes the average firing rate over [t0, tf] for each trial by drawing Mc samples from the smoothing distribution and evaluating the conditional intensity.

Parameters:

  • xK (array, shape (numBasis, K)) – Smoothed state estimates (basis coefficients × trials).

  • Wku (array, shape (numBasis, numBasis, K, K) or compatible) – Smoothed state covariance.

  • dN (array, shape (C, N)) – Observation (spike indicator) matrix.

  • t0 (float) – Time window over which to compute the average rate.

  • tf (float) – Time window over which to compute the average rate.

  • fitType (str) – 'poisson' or 'binomial'.

  • delta (float) – Time-step size in seconds.

  • gamma (array or None) – History-effect coefficients.

  • windowTimes (array or None) – History window boundaries.

  • Mc (int) – Number of Monte Carlo draws.

  • alphaVal (float) – Significance level for CIs (one-sided).

Returns:

  • spikeRateSig (Covariate) – Mean spike rate per trial with attached ConfidenceInterval.

  • ProbMat (ndarray, shape (K, K)) – ProbMat[k, m] = P(rate_m > rate_k) estimated from MC draws.

  • sigMat (ndarray, shape (K, K)) – Binary significance matrix at level 1 - alphaVal.

static computeSpikeRateDiffCIs(xK, Wku, dN, time1, time2, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, Mc=500, alphaVal=0.05)[source]

Monte Carlo CIs for the difference in spike rates between two time windows.

Computes the difference of average firing rates rate(time1) - rate(time2) for each trial by drawing Mc samples from the smoothing distribution.

Parameters:

  • xK (array, shape (numBasis, K)) – Smoothed state estimates (basis coefficients × trials).

  • Wku (array, shape (numBasis, numBasis, K, K) or compatible) – Smoothed state covariance.

  • dN (array, shape (C, N)) – Observation (spike indicator) matrix.

  • time1 (array-like, length 2) – [t0_1, tf_1] — first time window.

  • time2 (array-like, length 2) – [t0_2, tf_2] — second time window.

  • fitType (str) – 'poisson' or 'binomial'.

  • delta (float) – Time-step size in seconds.

  • gamma (array or None) – History-effect coefficients.

  • windowTimes (array or None) – History window boundaries.

  • Mc (int) – Number of Monte Carlo draws.

  • alphaVal (float) – Significance level for CIs (one-sided).

Returns:

  • spikeRateSig (Covariate) – Mean spike-rate difference per trial with attached CI.

  • ProbMat (ndarray, shape (K, K)) – ProbMat[k, m] = P(diff_m > diff_k) from MC draws.

  • sigMat (ndarray, shape (K, K)) – Binary significance matrix at level 1 - alphaVal.

static KF_EMCreateConstraints(EstimateA=1, AhatDiag=0, QhatDiag=1, QhatIsotropic=0, RhatDiag=1, RhatIsotropic=0, Estimatex0=1, EstimatePx0=1, Px0Isotropic=0, mcIter=1000, EnableIkeda=0)[source]

Return a dict of EM constraint flags for KF_EM().

Parameters:

  • EstimateA (int) – Whether to estimate the state transition matrix A.

  • AhatDiag (int) – Constrain A to be diagonal.

  • QhatDiag (int) – Constrain Q to be diagonal.

  • QhatIsotropic (int) – Constrain Q to be isotropic (scalar times identity). Only active when QhatDiag is also true.

  • RhatDiag (int) – Constrain R to be diagonal.

  • RhatIsotropic (int) – Constrain R to be isotropic. Only active when RhatDiag is also true.

  • Estimatex0 (int) – Whether to estimate the initial state x0.

  • EstimatePx0 (int) – Whether to estimate the initial covariance Px0.

  • Px0Isotropic (int) – Constrain Px0 to be isotropic. Only active when EstimatePx0 is true.

  • mcIter (int) – Number of Monte Carlo iterations for standard-error estimation via the observed information matrix.

  • EnableIkeda (int) – Enable Ikeda acceleration in the EM loop.

Returns:

Constraint dictionary consumed by KF_EM(), KF_MStep(), and KF_ComputeParamStandardErrors().

Return type:

dict

Notes

Docs only; mirrors nSTAT PR #138 (fix/kf-em @ 584ec96). QhatDiag=1 / RhatDiag=1 impose diagonal structure on Qhat / Rhat in the frame KF_EM internally whitens by the STARTING covariances Q0 / R0 (Tq = inv(chol(Q0, 'lower')), Tr = inv(chol(R0, 'lower')); see the G1 whitening fix), not in the caller’s own state / observation coordinates. For a diagonal Q0 / R0 the two frames coincide, so this is invisible. For a non-diagonal Q0 / R0, “Qhat/Rhat is diagonal” means Tq @ Qhat @ Tq.T / Tr @ Rhat @ Tr.T is diagonal; the returned Qhat / Rhat (mapped back to the caller’s coordinates) generally is not.

static KF_ComputeParamStandardErrors(y, xKFinal, WKFinal, Ahat, Qhat, Chat, Rhat, alphahat, x0hat, Px0hat, ExpectationSumsFinal, KFEM_Constraints=None)[source]

Compute standard errors via the observed information matrix.

Uses the complete information matrix and a Monte Carlo estimate of the missing information matrix, following McLachlan and Krishnan Eq. 4.7: Io(theta; y) = Ic(theta; y) - Im(theta; y).

Parameters:

  • y ((Dy, K) observations)

  • xKFinal ((Dx, K) smoothed states)

  • WKFinal ((Dx, Dx, K) smoothed covariances)

  • Ahat (estimated model matrices)

  • Qhat (estimated model matrices)

  • Chat (estimated model matrices)

  • Rhat (estimated model matrices)

  • alphahat ((Dy, 1) observation offset)

  • x0hat ((Dx,) initial state)

  • Px0hat ((Dx, Dx) initial covariance)

  • ExpectationSumsFinal (dict from KF_EStep())

  • KFEM_Constraints (dict from KF_EMCreateConstraints())

Returns:

  • SE (dict of standard-error matrices/vectors for each parameter)

  • Pvals (dict of p-value matrices/vectors for each parameter)

static KF_EM(y, Ahat0, Qhat0, Chat0, Rhat0, alphahat0, x0=None, Px0=None, KFEM_Constraints=None)[source]

Kalman Filter EM algorithm with Cholesky-scaled system.

Estimates the parameters of a linear-Gaussian state-space model:

x_{k+1} = A x_k + v_k,  v ~ N(0, Q)
y_k     = C x_k + alpha + w_k,  w ~ N(0, R)

using the Expectation-Maximisation algorithm (E-step: KF + RTS smoother, M-step: closed-form updates). Optionally applies Ikeda acceleration.

Parameters:

  • y ((Dy, K) observation matrix (each column is one time step))

  • Ahat0 ((Dx, Dx) initial state transition)

  • Qhat0 ((Dx, Dx) initial process noise covariance)

  • Chat0 ((Dy, Dx) initial observation matrix)

  • Rhat0 ((Dy, Dy) initial observation noise covariance)

  • alphahat0 ((Dy, 1) initial observation offset)

  • x0 ((Dx,) initial state (default zeros))

  • Px0 ((Dx, Dx) initial state covariance (default 1e-10 * I))

  • KFEM_Constraints (dict from KF_EMCreateConstraints())

Returns:

  • xKFinal ((Dx, K) smoothed state estimates)

  • WKFinal ((Dx, Dx, K) smoothed covariances)

  • Ahat, Qhat, Chat, Rhat (estimated model matrices)

  • alphahat ((Dy, 1) estimated observation offset)

  • x0hat ((Dx,) estimated initial state)

  • Px0hat ((Dx, Dx) estimated initial covariance)

  • IC (dict of information criteria (AIC, AICc, BIC, llobs, llcomp))

  • SE (dict of standard errors (or empty dict if not computed))

  • Pvals (dict of p-values (or empty dict if not computed))

  • nIter (int — number of EM iterations)

Notes

Docs only; mirrors nSTAT PR #138 (fix/kf-em @ 584ec96 / dd7c7d5). EM stops at the FIRST decrease of the log-likelihood, not at convergence; see PP_EM()’s identical note for the general rule. Unlike PP_EM / PPLFP_EM (whose Monte Carlo M-step makes the stopping iteration random run to run), KF_EM’s E-step and M-step are both closed-form and deterministic, so its stopping iteration is NOT random for a fixed problem – though it still depends on the starting parameters Ahat0/Qhat0/…/Px0.

static KF_EStep(A, Q, C, R, y, alpha, x0, Px0)[source]

E-step for the Kalman Filter EM algorithm.

Runs the forward Kalman filter followed by the backward RTS smoother and computes sufficient statistics (expectation sums) for the M-step.

Parameters:

  • A ((Dx, Dx) state transition matrix)

  • Q ((Dx, Dx) process noise covariance)

  • C ((Dy, Dx) observation matrix)

  • R ((Dy, Dy) observation noise covariance)

  • y ((Dy, K) observations (state-major))

  • alpha ((Dy, 1) or (Dy,) observation offset)

  • x0 ((Dx,) initial state)

  • Px0 ((Dx, Dx) initial state covariance)

Returns:

  • x_K ((Dx, K) smoothed states)

  • W_K ((Dx, Dx, K) smoothed covariances)

  • logll (float — complete-data log-likelihood)

  • ExpectationSums (dict of sufficient statistics)

static KF_MStep(y, x_K, x0, Px0, ExpectationSums, KFEM_Constraints=None)[source]

M-step for the Kalman Filter EM algorithm.

Updates all state-space model parameters given the sufficient statistics from KF_EStep().

Parameters:

  • y ((Dy, K) observations)

  • x_K ((Dx, K) smoothed states from E-step)

  • x0 ((Dx,) current initial state estimate)

  • Px0 ((Dx, Dx) current initial covariance estimate)

  • ExpectationSums (dict from KF_EStep())

  • KFEM_Constraints (dict from KF_EMCreateConstraints(), or None)

Return type:

Ahat, Qhat, Chat, Rhat, alphahat, x0hat, Px0hat

static PPLFP_fixedIntervalSmoother(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_fixedIntervalSmoother instead.

static PPLFP_DecodeLinear(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_DecodeLinear instead.

static PPLFP_Decode_predict(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_Decode_predict instead.

static PPLFP_Decode_update(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_Decode_update instead.

static PPLFP_EMCreateConstraints(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_EMCreateConstraints instead.

static PPLFP_ComputeParamStandardErrors(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_ComputeParamStandardErrors instead.

static PPLFP_EM(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_EM instead.

static PPLFP_EStep(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_EStep instead.

static PPLFP_MStep(*args, **kwargs)[source]

[DEPRECATED] Use nstat.decoding.PPLFP.PPLFP_MStep instead.

static PP_EMCreateConstraints(EstimateA=1, AhatDiag=0, QhatDiag=1, QhatIsotropic=0, Estimatex0=0, EstimatePx0=0, Px0Isotropic=0, mcIter=1000, EnableIkeda=0)[source]

Build a constraints dict for PP_EM.

Matlab: PointProcessEM.PP_EMCreateConstraints. Defaults follow the repaired MATLAB (fix/pp-em round 2): EstimateA=1, AhatDiag=0, QhatDiag=1, QhatIsotropic=0, Estimatex0=0, EstimatePx0=0, Px0Isotropic=0, mcIter=1000, EnableIkeda=0. x0 and Px0 are NOT estimated by default (they were 1 and 1): the Px0 M-step, Px0hat = (x0hat - x0)(x0hat - x0)' .* I, is a single-sample estimate that collapses to ~0 after one iteration, which drives -1/2 log det(Px0) and hence the E-step log-likelihood to +Inf and stops EM after ~2 iterations (MATLAB report C2). Pass Estimatex0=1, EstimatePx0=1 explicitly for the old behaviour.

Parameters:

  • EstimateA (int) – Whether to estimate the state transition matrix A.

  • AhatDiag (int) – Constrain A to be diagonal.

  • QhatDiag (int) – Constrain Q to be diagonal.

  • QhatIsotropic (int) – Constrain Q to be isotropic (scalar * I).

  • Estimatex0 (int) – Whether to estimate the initial state x0 (default 0; see above).

  • EstimatePx0 (int) – Whether to estimate the initial state covariance Px0 (default 0).

  • Px0Isotropic (int) – Constrain Px0 to be isotropic.

  • mcIter (int) – Number of Monte Carlo iterations for standard error estimation.

  • EnableIkeda (int) – Enable Ikeda acceleration. Not ported: PP_EM raises NotImplementedError when it is 1 (MATLAB runs an Ikeda step, history-free models only).

Returns:

Constraints dictionary with all fields.

Return type:

dict

Notes

Docs only; mirrors nSTAT PR #138 (fix/kf-em @ 584ec96). QhatDiag=1 imposes diagonal structure on Qhat in the frame PP_EM internally whitens by the STARTING covariance Q0 (Tq = inv(chol(Q0, 'lower')); see G1), not in the caller’s own state coordinates. For a diagonal Q0 the two frames coincide, so this is invisible. For a non-diagonal Q0, “Qhat is diagonal” means Tq @ Qhat @ Tq.T is diagonal; the returned Qhat (mapped back to the caller’s coordinates) generally is not.

static PP_ComputeParamStandardErrors(dN, xKFinal, WKFinal, Ahat, Qhat, x0hat, Px0hat, ExpectationSumsFinal, fitType, muhat, betahat, gammahat, windowTimes, HkAll, PPEM_Constraints=None)[source]

Compute standard errors via the observed information matrix.

Uses a Monte-Carlo approximation of the missing information matrix (McLachlan & Krishnan, Eq. 4.7), with mcIter draws from NumPy’s global stream (seed with np.random.seed or nstat.extras.matlab_rng.seeded_global_rng()). As in MATLAB (nSTAT PR #135): the constraints are honoured, the intensities are unclipped, the inverse observed information is projected with MATLAB’s nearestSPD and the p-values are MATLAB ztest’s 2*normcdf(-|z|). An exactly singular observed information (MATLAB’s SE pass never returns there) is a Python-only extension: the SEs come from its pseudo-inverse, every parameter with a component in its null space – not identifiable from the data, typically the coefficient of a separated history window at the exp() underflow – gets SE and p-value NaN, and a RuntimeWarning names those parameters.

Parameters:

  • dN ((C, N) spike observations)

  • xKFinal ((dx, N) smoothed states)

  • WKFinal ((dx, dx, N) smoothed state covariances)

  • Ahat ((dx, dx) estimated state transition)

  • Qhat ((dx, dx) estimated state noise covariance)

  • x0hat ((dx,) estimated initial state)

  • Px0hat ((dx, dx) estimated initial state covariance)

  • ExpectationSumsFinal (dict with sufficient statistics)

  • fitType ('poisson' or 'binomial')

  • muhat ((C,) estimated baseline rates)

  • betahat ((dx, C) estimated stimulus coefficients)

  • gammahat ((nW, C) or scalar estimated history coefficients)

  • windowTimes (history window boundaries or None)

  • HkAll ((N, nW, C) history design tensor)

  • PPEM_Constraints (dict from PP_EMCreateConstraints)

Returns:

  • SE (dict of standard errors for each parameter group)

  • Pvals (dict of p-values for each parameter group)

  • nTerms (int, total number of estimated parameters)

static PP_EM(dN, Ahat0, Qhat0, mu, beta, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, x0=None, Px0=None, PPEM_Constraints=None, MstepMethod='NewtonRaphson')[source]

Full Point-Process state-space EM algorithm.

Estimates state-space model parameters (A, Q, mu, beta, gamma) via EM for point-process observations. Unlike PPSS_EM, this operates on raw spike observations with explicit beta/mu/gamma parameters (no basis functions).

Matlab: PointProcessEM.PP_EM (repaired, fix/pp-em). Defaults: fitType='poisson', delta=0.001, no history (gamma empty or 0 with windowTimes empty), x0 = 0, Px0 = 1e-9 I, PPEM_Constraints = PP_EMCreateConstraints() (x0 / Px0 not estimated: that estimator collapses Px0 and stops EM after ~2 iterations) and MstepMethod='NewtonRaphson' – the Monte-Carlo Newton-Raphson M-step maximises the expected complete-data log-likelihood over the smoothed state posterior. (The alternative ‘GLM’ M-step is a plug-in regression on the smoothed means that ignores W_K, inflates beta and drifts; it is still available by passing it explicitly. Any other MstepMethod raises ValueError.) With gamma given and windowTimes empty there is one history window per coefficient, windowTimes = 0:delta:size(gamma,1)*delta, and a nonzero shared numWindows x 1 gamma is expanded to every cell. History is built on the delta grid. EM stops before the M-step if the E-step log-likelihood is not finite and returns the best finite iterate. EnableIkeda = 1 (MATLAB’s Ikeda acceleration step) is not ported and raises NotImplementedError; a Qhat0 that is not positive definite raises LinAlgError (MATLAB’s chol errors). The Monte Carlo of the M-steps and the SE pass draws from NumPy’s global stream, so np.random.seed(s) (or nstat.extras.matlab_rng.seeded_global_rng()) makes a fit reproducible; MATLAB’s randn stream is not reproduced (the end-to-end gold em_drivers.mat is compared within the measured Monte Carlo spread).

SE / Pvals are always computed (MATLAB computes them only when more than 10 outputs are requested). The Monte Carlo SE pass dominates the run time at the default mcIter = 1000: measured 4.5 s of 6.2 s for N = 800 bins, C = 4 cells, 3 history windows (0.5 s of 2.1 s at mcIter = 100); pass a smaller mcIter in PPEM_Constraints to reduce it.

Note (docs only; mirrors nSTAT PR #138, fix/kf-em @ 584ec96): EM stops at the FIRST decrease of the log-likelihood (dLikelihood < 0), not at convergence. Under the Newton-Raphson M-step’s Monte Carlo expectation, successive iterations’ log-likelihoods are themselves noisy, so the iteration at which that first decrease happens – and hence nIter, and which iterate xKFinal/Ahat/…/IC come from – is random (the same problem, run repeatedly, was observed to stop anywhere from iteration 6 to 11 on one problem). This is a property of the stopping rule, not a bug; a smaller mcIter increases the noise and the variability. KF_EM()’s closed-form E/M steps make its stopping iteration deterministic for a fixed problem, though it still depends on the starting parameters; see its docstring.

Parameters:

  • dN ((C, N) binary spike observations)

  • Ahat0 ((dx, dx) initial state transition matrix)

  • Qhat0 ((dx, dx) initial state noise covariance)

  • mu ((C,) initial baseline log-rates)

  • beta ((dx, C) initial stimulus coefficients)

  • fitType ('poisson' or 'binomial')

  • delta (float, time bin width)

  • gamma ((nW, C) or scalar, initial history coefficients)

  • windowTimes (history window boundaries)

  • x0 ((dx,) initial state (default zeros))

  • Px0 ((dx, dx) initial state covariance)

  • PPEM_Constraints (dict from PP_EMCreateConstraints)

  • MstepMethod ('NewtonRaphson' or 'GLM' (passed to PP_MStep with) – delta)

Returns:

  • xKFinal, WKFinal, Ahat, Qhat, muhat, betahat, gammahat,

  • x0hat, Px0hat, IC, SE, Pvals, nIter

static PP_EStep(A, Q, dN, mu, beta, fitType, gamma, HkAll, x0, Px0)[source]

E-step for PP EM: forward filter + RTS smoother + cross-covariance.

Parameters:

  • A ((dx, dx) state transition matrix)

  • Q ((dx, dx) state noise covariance)

  • dN ((C, N) binary spike observations)

  • mu ((C,) baseline log-rate)

  • beta ((dx, C) stimulus coefficients)

  • fitType ('poisson' or 'binomial')

  • gamma ((nW, C) or scalar history coefficients)

  • HkAll ((N, nW, C) history design tensor)

  • x0 ((dx,) initial state)

  • Px0 ((dx, dx) initial state covariance)

Returns:

  • x_K ((dx, N) smoothed states)

  • W_K ((dx, dx, N) smoothed covariances)

  • logll (float, log-likelihood)

  • ExpectationSums (dict of sufficient statistics)

Notes

As in MATLAB the intensity is the unclipped exp(terms) (binomial exp./(1+exp)) and log(det(Q)) / log(det(Px0)) are unfloored: a non-positive determinant gives a non-finite logll, on which PP_EM stops before the M-step.

static PP_MStep(dN, x_K, W_K, x0, Px0, ExpectationSums, fitType, muhat, betahat, gammahat, windowTimes, HkAll, PPEM_Constraints=None, MstepMethod='NewtonRaphson', delta=0.001, *, _warn_glm=True)[source]

M-step for PP EM: update all model parameters.

Matlab: PointProcessEM.PP_MStep (repaired, fix/pp-em @ 8dbd0e4). A, Q, x0 and Px0 have closed-form updates; (mu, beta, gamma) are updated by the Monte-Carlo Newton-Raphson step (MstepMethod = 'NewtonRaphson', the default) or by a GLM fit on the smoothed means ('GLM'). Any other MstepMethod raises ValueError (MATLAB runs Newton-Raphson for it).

As in MATLAB (nSTAT PR #135, aa88a2b): the closed-form updates have no ridges or eigenvalue floors (a single-sample Px0hat can collapse to 0); each Newton loop runs at most 99 steps on the unclipped exp(terms) (binomial exp./(1+exp)), each step is MATLAB’s H\g (g/H for one coefficient, MATLAB’s LU dividing by the pivots for several: a denormal pivot, as on a separated history window near the exp() underflow, gives MATLAB’s finite step); as in MATLAB, a step with a NaN keeps the previous value and an infinite one is taken (an exactly singular Hessian gives either, e.g. -[1 1; 1 1]\[1; 2] = [-Inf; Inf]). The McExp = 50 draws come from NumPy’s global stream (seed with np.random.seed or nstat.extras.matlab_rng.seeded_global_rng()).

Parameters:

  • dN ((C, N) spike observations)

  • x_K ((dx, N) smoothed states)

  • W_K ((dx, dx, N) smoothed covariances)

  • x0 ((dx,) current initial state estimate)

  • Px0 ((dx, dx) current initial covariance estimate)

  • ExpectationSums (dict from E-step)

  • fitType ('poisson' or 'binomial')

  • muhat ((C,) current baseline rates)

  • betahat ((dx, C) current stimulus coefficients)

  • gammahat (scalar or (nW, C) current history coefficients)

  • windowTimes (history window boundaries or None)

  • HkAll ((N, nW, C) history tensor)

  • PPEM_Constraints (dict from PP_EMCreateConstraints (default) – PP_EMCreateConstraints(): x0 / Px0 not estimated)

  • MstepMethod ('NewtonRaphson' (default, as in the repaired MATLAB) or) – ‘GLM’. ‘GLM’ regresses each cell’s spikes on the smoothed means x_K (a plug-in fit that ignores W_K; see _em_glm_mstep()), as MATLAB’s GLM M-step does.

  • delta (float, default 0.001) – Seconds per bin (MATLAB’s optional 15th input): the GLM M-step’s time base. PP_EM passes its delta.

Return type:

Ahat, Qhat, muhat_new, betahat_new, gammahat_new, x0hat, Px0hat

static mPPCO_fixedIntervalSmoother(A, Q, C, R, y, alpha, dN, lags, mu, beta, fitType, delta=0.001, gamma=None, windowTimes=None, x0=None, Px0=None, HkAll=None)[source]

[DEPRECATED] Alias of PPLFP_fixedIntervalSmoother() (fixed-interval lag smoother).

Matlab: DecodingAlgorithms.mPPCO_fixedIntervalSmoother is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_fixedIntervalSmoother (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_fixedIntervalSmoother(), so it returns exactly what PPLFP_fixedIntervalSmoother returns: (x_pLag, W_pLag, x_uLag, W_uLag).

Use PPLFP_fixedIntervalSmoother() (or nstat.decoding.PPLFP.PPLFP_fixedIntervalSmoother) instead.

static mPPCODecodeLinear(A, Q, C, R, y, alpha, dN, mu, beta, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, x0=None, Px0=None, HkAll=None)[source]

[DEPRECATED] Alias of PPLFP_DecodeLinear() (the PPLFP forward filter).

Matlab: DecodingAlgorithms.mPPCODecodeLinear is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_DecodeLinear (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_DecodeLinear(), so it returns exactly what PPLFP_DecodeLinear returns: (x_p, W_p, x_u, W_u).

The legacy Python signature defaults fitType to 'poisson', which PPLFP_DecodeLinear (and MATLAB) do not; the signature is frozen, so that default is forwarded as-is. delta=0.001 equals PPLFP_DecodeLinear’s None default, which it resolves to 0.001.

Use PPLFP_DecodeLinear() (or nstat.decoding.PPLFP.PPLFP_DecodeLinear) instead.

static mPPCODecode_predict(x_u, W_u, A, Q)[source]

[DEPRECATED] Alias of PPLFP_Decode_predict() (the PPLFP predict step).

Matlab: DecodingAlgorithms.mPPCODecode_predict is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_Decode_predict (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_Decode_predict(), so it returns exactly what PPLFP_Decode_predict returns: (x_p, W_p) with x_p = A @ x_u (an (ns, 1) column x_u now gives an (ns, 1) x_p, as in MATLAB, where the old body flattened it).

Use PPLFP_Decode_predict() (or nstat.decoding.PPLFP.PPLFP_Decode_predict) instead.

static mPPCODecode_update(x_p, W_p, C, R, y, alpha, dN, mu, beta, fitType='poisson', gamma=None, HkAll=None, time_index=0, WuConv=None)[source]

[DEPRECATED] Alias of PPLFP_Decode_update() (the PPLFP measurement update).

Matlab: DecodingAlgorithms.mPPCODecode_update (lines 4855-4944) is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_Decode_update (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards to PPLFP_Decode_update(), so it returns exactly what PPLFP_Decode_update returns: (x_u, W_u, lambdaDeltaMat).

Unlike the other mPPCO_* aliases, this one’s documented MATLAB contract takes HkAll permuted as (numWindows, numCells, N) (time on the 3rd axis) – the layout this alias’s own standalone body used before this fix, and the layout callers of this frozen API signature still pass. PPLFP_Decode_update takes the canonical (N, numWindows, numCells) layout, so HkAll is transposed (np.transpose(HkAll, (2, 0, 1))) before forwarding; every other argument forwards unchanged.

Use PPLFP_Decode_update() (or nstat.decoding.PPLFP.PPLFP_Decode_update, with the canonical (N, numWindows, numCells) HkAll) instead.

static mPPCO_EMCreateConstraints(EstimateA=1, AhatDiag=0, QhatDiag=1, QhatIsotropic=0, RhatDiag=1, RhatIsotropic=0, Estimatex0=0, EstimatePx0=0, Px0Isotropic=0, mcIter=1000, EnableIkeda=0)[source]

[DEPRECATED] Alias of PPLFP_EMCreateConstraints() (EM constraints builder).

Matlab: DecodingAlgorithms.mPPCO_EMCreateConstraints is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_EMCreateConstraints (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_EMCreateConstraints(), so it returns exactly what PPLFP_EMCreateConstraints returns: the constraints dict.

Its defaults are PPLFP_EMCreateConstraints’s (MATLAB forwards varargin, so the alias inherits them): Estimatex0=0 and EstimatePx0=0 since the repaired MATLAB (B8).

Use PPLFP_EMCreateConstraints() (or nstat.decoding.PPLFP.PPLFP_EMCreateConstraints) instead.

static mPPCO_ComputeParamStandardErrors(y, dN, xKFinal, WKFinal, Ahat, Qhat, Chat, Rhat, alphahat, x0hat, Px0hat, ExpectationSumsFinal, fitType, muhat, betahat, gammahat, windowTimes, HkAll, mPPCOEM_Constraints=None)[source]

[DEPRECATED] Alias of PPLFP_ComputeParamStandardErrors() (parameter standard errors).

Matlab: DecodingAlgorithms.mPPCO_ComputeParamStandardErrors is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_ComputeParamStandardErrors (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_ComputeParamStandardErrors(), so it returns exactly what PPLFP_ComputeParamStandardErrors returns: (SE, Pvals, nTerms). mPPCOEM_Constraints is the legacy name of PPLFP_EM_Constraints; it is passed positionally.

Use PPLFP_ComputeParamStandardErrors() (or nstat.decoding.PPLFP.PPLFP_ComputeParamStandardErrors) instead.

static mPPCO_EM(y, dN, Ahat0, Qhat0, Chat0, Rhat0, alphahat0, mu, beta, fitType='poisson', delta=0.001, gamma=None, windowTimes=None, x0=None, Px0=None, mPPCOEM_Constraints=None, MstepMethod='NewtonRaphson')[source]

[DEPRECATED] Alias of PPLFP_EM() (the full EM driver).

Matlab: DecodingAlgorithms.mPPCO_EM is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_EM (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_EM(), so it returns exactly what PPLFP_EM returns: (xKFinal, WKFinal, Ahat, Qhat, Chat, Rhat, alphahat, muhat, betahat, gammahat, x0hat, Px0hat, IC, SE, Pvals). mPPCOEM_Constraints is the legacy name of PPLFP_EM_Constraints; it is passed positionally. MstepMethod defaults to PPLFP_EM’s 'NewtonRaphson' (the repaired MATLAB’s, B8, which the MATLAB alias inherits through varargin).

Use PPLFP_EM() (or nstat.decoding.PPLFP.PPLFP_EM) instead.

static mPPCO_EStep(A, Q, C, R, y, alpha, dN, mu, beta, fitType='poisson', delta=0.001, gamma=None, HkAll=None, x0=None, Px0=None)[source]

[DEPRECATED] Alias of PPLFP_EStep() (the EM expectation step).

Matlab: DecodingAlgorithms.mPPCO_EStep is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_EStep (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_EStep(), so it returns exactly what PPLFP_EStep returns: (x_K, W_K, logll, ExpectationSums).

Use PPLFP_EStep() (or nstat.decoding.PPLFP.PPLFP_EStep) instead.

static mPPCO_MStep(dN, y, x_K, W_K, x0, Px0, ExpectationSums, fitType='poisson', muhat=None, betahat=None, gammahat=None, windowTimes=None, HkAll=None, mPPCOEM_Constraints=None, MstepMethod='NewtonRaphson', delta=0.001)[source]

[DEPRECATED] Alias of PPLFP_MStep() (the EM maximisation step).

Matlab: DecodingAlgorithms.mPPCO_MStep is a deprecation shim that warns nSTAT:deprecated:mPPCO and forwards varargin{:} to DecodingAlgorithms.PPLFP_MStep (the mPPCO family was renamed PPLFP). This mirror emits a DeprecationWarning with MATLAB’s message and forwards every argument positionally to PPLFP_MStep(), so it returns exactly what PPLFP_MStep returns: (Ahat, Qhat, Chat, Rhat, alphahat, muhat_new, betahat_new, gammahat_new, x0hat, Px0hat). mPPCOEM_Constraints is the legacy name of PPLFP_EM_Constraints; it is passed positionally.

Unlike the other aliases, the legacy Python signature gives fitType / muhat / betahat / gammahat / windowTimes / HkAll defaults that PPLFP_MStep (and MATLAB) do not have; the signature is frozen, so those defaults are forwarded as-is. Supply all six: omitting them is unsupported (MATLAB raises, here PPLFP_MStep receives None). MstepMethod defaults to 'NewtonRaphson', PPLFP_MStep’s default since the repaired MATLAB (B8), which the MATLAB alias inherits through varargin. The trailing delta (default 0.001) is PPLFP_MStep’s optional 16th input of the repaired MATLAB (R4c), which the MATLAB alias also forwards through varargin.

Use PPLFP_MStep() (or nstat.decoding.PPLFP.PPLFP_MStep) instead.