neurotools.stats.glm module

Tutorial in Poisson generalized linear point-process models for neural spike train analysis.

Depends on

neurotools.stats.glm.GLMPenaltyPoisson(X, Y)[source]

Generates objective, gradient, and hessian functions for the Poisson GLM for design matrix X and spike observations Y.

Parameters:
  • X (array-like) – N observation x D features design matrix

  • Y (array-like) – N x 1 point-process counts

Returns:

  • objective (function) – objective function for scipy.optimize.minimize

  • gradient (function) – gradient (jacobian) function for scipy.optimize.minimize

  • hessian (function) – hessian function for scipy.optimize.minimize

neurotools.stats.glm.sexp(x)[source]

Exponential function avoiding overflow and underflow, for evaluating the logistic sigmoid function

neurotools.stats.glm.sigmoid(x)[source]

Logistic sigmoid function

neurotools.stats.glm.GLMPenaltyBernoulli(X, Y)[source]

Generates objective, gradient, and hessian functions for the Bernoulli GLM for design matrix X and spike observations Y.

Parameters:
  • X (array-like) – N observation x D features design matrix

  • Y (array-like) – N x 1 point-process counts

Returns:

  • objective (function) – objective function for scipy.optimize.minimize

  • gradient (function) – gradient (jacobian) function for scipy.optimize.minimize

  • hessian (function) – hessian function for scipy.optimize.minimize

neurotools.stats.glm.GLMPenaltyL2(X, Y, penalties=None)[source]

Generates objective, gradient, and hessian functions for the penalized L2 regularized poisson GLM for design matrix X and spike observations Y.

Parameters:
  • X (array-like) – N observation x D features design matrix

  • Y (array-like) – N x 1 point-process counts

  • penalties – len D-1 list of L2 penalty weights (don’t penalize mean)

Returns:

  • objective (function) – objective function for scipy.optimize.minimize

  • gradient (function) – gradient (jacobian) function for scipy.optimize.minimize

  • hessian (function) – hessian function for scipy.optimize.minimize

neurotools.stats.glm.ppglmfit(X, Y, verbose=False)[source]

The GLM solver in statsmodels is very general. It accepts any link function and expects that, if you want a constant term in your model, that you have already manually added a column of ones to your design matrix. This wrapper simplifies using GLM to fit the common case of a Poisson point-process model, where the constant term has not been explicitly added to the design matrix

Parameters:
  • X (N_observations x N_features design matrix.)

  • Y (Binary point process observations)

Returns:

μ, B

Return type:

the offset and parameter estimates for the GLM model.

neurotools.stats.glm.fitGLM(X, Y, L2Penalty=0.0, mu0=None, B0=None, **kwargs)[source]

Fit the model using gradient descent with hessian

Parameters:
  • X (matrix) – design matrix

  • Y (vector) – binary spike observations

  • L2Penalty (scalar) – optional L2 penalty on features, defaults to 0

  • mu0 (scalar or None) – Initial guess for mean offset. Optional, defaults to 0 if None.

  • B0 (vector or None) – Initial guess for parameter weights. Optional, defaults to zeroes if None

  • kwargs – Keyword arguments to forward to scipy.optimize.minimize

Returns:

  • mu (mean offset parameter)

  • B (feature weights)

neurotools.stats.glm.crossvalidatedAUC(X, Y, NXVAL=4)[source]

Crossvalidated area under the ROC curve calculation. This routine uses the non-regularized GLMPenaltyL2 to fit a GLM point-process model and test accuracy under K-fold crossvalidation.

Parameters:
  • X (np.array) – Covariate matrix Nsamples x Nfeatures

  • Y (np.array) – Binary point-process observations, 1D array length Nsamples

  • NXVAL (positive int) – Defaults to 4. Number of cross-validation blocks to use

Returns:

Area under the ROC curve, cross-validated, for non-regularized GLM point process model fit

Return type:

float

neurotools.stats.glm.gradientglmfit(X, Y, L2Penalty=0.0)[source]

mu_hat, B_hat = gradientglmfit(X,Y,L2Penalty=0.0)

Fit Poisson GLM using gradient descent with hessian

Parameters:
  • X (np.array) – Covariate matrix Nsamples x Nfeatures

  • Y (np.array) – Binary point-process observations, 1D array length Nsamples

neurotools.stats.glm.cosine_kernel(x)[source]

Raised cosine basis kernel, normalized such that it integrates to 1 centered at zero. Time is rescaled so that the kernel spans from -2 to 2

Parameters:

x (vector)

Returns:

$tfrac 1 4 + tfrac 1 2 cos(x)$ if $xin[-pi,pi]$, otherwise 0.

Return type:

vector

neurotools.stats.glm.log_cosine_basis(N=range(1, 6), t=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]), base=2, offset=1, normalize=True)[source]

Generate overlapping log-cosine basis elements

Parameters:
  • N (array) – Array of wave quarter-phases

  • t (array) – times

  • base (scalar) – exponent base

  • offset (scalar) – leave this set to 1 (default)

Returns:

B – Basis with n_elements x n_times shape

Return type:

array

neurotools.stats.glm.make_cosine_basis(N, L, min_interval, normalize=True)[source]

Build N logarightmically spaced cosine basis functions spanning L samples, with a peak resolution of min_interval

# Solve for a time basis with these constraints # t[0] = 0 # t[min_interval] = 1 # log(L)/log(b) = n_basis+1 # log(b) = log(L)/(n_basis+1) # b = exp(log(L)/(n_basis+1))

Parameters:
  • N (int) – Number of basis functions

  • L (int) – Number of time-bins

  • min_interval (scalar) – Number of bins between the two shortes basis elements. That is, minimum time separation between basis functions.

Returns:

B – Basis with n_elements x n_times shape

Return type:

array

neurotools.stats.glm.numeric_grad(obj, p, delta=1.1920929e-07)[source]

Numerically estimate the gradient of an objective function

neurotools.stats.glm.numeric_hess(jac, p, delta=1.1920929e-07)[source]

Numerically estimate the hessian given a gradient function