use matRegress instead of glmfit

it is often faster, and we have full control over the details
Note that matRegress is a subtree of neuroGLM (careful not to mix histories!)

tutorial.m

@@ -78,14 +78,33 @@
 colIndices = buildGLM.getDesignMatrixColIndices(dspec, 'LFP');
 dm = buildGLM.zscoreDesignMatrix(dm, [colIndices{:}]);
 
-%% Maximum likelihood estimation using glmfit
-[w, dev, stats] = glmfit(dm.X, y, 'poisson', 'link', 'log');
-DC = w(1); w = w(2:end);
-DC_se = stats.se(1); stats.se = stats.se(2:end);
+dm = buildGLM.addBiasColumn(dm); % DO NOT ADD THE BIAS TERM IF USING GLMFIT
+
+%% Least squares for initialization
+tic
+wInit = dm.X \ y;
+toc
+
+%% Use matRegress for Poisson regression
+% it requires `fminunc` from MATLAB's optimization toolbox
+addpath('matRegress')
+
+fnlin = @nlfuns.exp; % inverse link function (a.k.a. nonlinearity)
+lfunc = @(w)(glms.neglog.poisson(w, dm.X, y, fnlin)); % cost/loss function
+
+opts = optimoptions(@fminunc, 'Algorithm', 'trust-region', ...
+    'GradObj', 'on', 'Hessian','on');
+
+[wml, nlogli, exitflag, ostruct, grad, hessian] = fminunc(lfunc, wInit, opts);
+wvar = diag(inv(hessian));
+
+%% Alternative maximum likelihood Poisson estimation using glmfit
+% [w, dev, stats] = glmfit(dm.X, y, 'poisson', 'link', 'log');
+% wvar = stats.se.^2;
 
 %% Visualize
-ws = buildGLM.combineWeights(dm, w);
-wvar = buildGLM.combineWeights(dm, stats.se.^2);
+ws = buildGLM.combineWeights(dm, wml);
+wvar = buildGLM.combineWeights(dm, wvar);
 
 fig = figure(2913); clf;
 nCovar = numel(dspec.covar);