Bias Column handling

+basisFactory/deltaStim.m

@@ -1,6 +1,7 @@
 function stim = deltaStim(bt, nT, v)
 % Returns a sparse vector with events at binned timings
 
+bt(bt>nT) = [];
 o = ones(numel(bt), 1);
 
 if nargin < 3

+basisFactory/makeNonlinearRaisedCos.m

@@ -1,4 +1,4 @@
-function bases = makeNonlinearRaisedCos(nBases, binSize, endPoints, nlOffset)
+function bases = makeNonlinearRaisedCos(nBases, binSize, endPoints, nlOffset, varargin)
 % Make nonlinearly stretched basis consisting of raised cosines.
 % Nonlinear stretching allows faster changes near the event.
 %
@@ -8,6 +8,10 @@
 %          (i.e. center) of the last raised cosine basis vectors
 %   nlOffset: [1] offset for nonlinear stretching of x axis:  y = log(t+nlOffset)
 %         (larger nlOffset -> more nearly linear stretching)
+
+%   Optional Arguments (entered as argument pairs)
+%       'Normalize' (default = false) normalizes basis vectors to sum to 1   
+%       'Orthogonalize' (default = false) orthogonalizes the basis so B'*B=1
 %
 %  Outputs:  iht = time lattice on which basis is defined
 %            ihbasis = basis itself
@@ -16,6 +20,11 @@
 %  Example call
 %  bases = basisFactory.makeNonlinearRaisedCos(10, 1, [0 500], 2);
 
+p = inputParser();
+p.addOptional('Normalize', false);
+p.addOptional('Orthogonalize', false);
+p.parse(varargin{:});
+
 % nonlinearity for stretching x axis (and its inverse)
 nlin = @(x)(log(x + 1e-20));
 invnl = @(x)(exp(x) - 1e-20);
@@ -33,12 +42,20 @@
 ihbasis = ff(repmat(nlin(iht + nlOffset), 1, nBases), repmat(ctrs, numel(iht), 1), db);
 ihctrs = invnl(ctrs);
 
+if p.Results.Normalize
+    ihbasis = bsxfun(@rdivide, ihbasis, sum(ihbasis));
+end
+
+if p.Results.Orthogonalize
+    ihbasis = orth(ihbasis);
+end
+
 bases.type = mfilename;
 bases.param.nBases = nBases;
 bases.param.binSize = binSize;
 bases.param.endPoints = endPoints;
 bases.param.nlOffset = nlOffset;
 bases.B = ihbasis;
 bases.edim = size(bases.B, 2);
-bases.tr = iht;
+bases.tr = iht/binSize;
 bases.centers = ihctrs;

+basisFactory/makeSmoothTemporalBasis.m

@@ -1,47 +1,64 @@
-function bases = makeSmoothTemporalBasis(shape, duration, nBases, binfun)
+function bases = makeSmoothTemporalBasis(shape, duration, nBases, binfun, varargin)
 %
 % Input
 %   shape: 'raised cosine' or 'boxcar'
 %   duration: the time that needs to be covered
 %   nBases: number of basis vectors to fill the duration
-%   binfun: 
+%   binfun:
+%
+%   Optional Arguments (entered as argument pairs)
+%       'Normalize' (default = false) normalizes basis vectors to sum to 1
+%       'Orthogonalize' (default = false) orthogonalizes the basis so B'*B=1
 %
 % Output
 %   BBstm: basis vectors
 
+p = inputParser();
+p.addOptional('Normalize', false);
+p.addOptional('Orthogonalize', false);
+p.parse(varargin{:});
+
 nkbins = binfun(duration); % number of bins for the basis functions
 
 ttb = repmat((1:nkbins)', 1, nBases); % time indices for basis
 
-if strcmpi(shape, 'raised cosine')
-    %   ^
-    %  / \
-    % /   \______
-    %      ^
-    %     / \
-    % ___/   \___
-    %         ^
-    %        / \
-    % ______/   \
-    % For raised cosine, the spacing between the centers must be 1/4 of the
-    % width of the cosine
-    dbcenter = nkbins / (3 + nBases); % spacing between bumps
-    width = 4 * dbcenter; % width of each bump
-    bcenters = 2 * dbcenter + dbcenter*(0:nBases-1);
-    % location of each bump centers
-    % bcenters = round(bcenters); % round to closest ms <-- BAD!!! NEVER DO THIS
-    bfun = @(x,period)((abs(x/period)<0.5).*(cos(x*2*pi/period)*.5+.5));
-    BBstm = bfun(ttb-repmat(bcenters,nkbins,1), width);
-elseif strcmpi(shape, 'boxcar')
-    width = nkbins / nBases;
-    BBstm = zeros(size(ttb));
-    bcenters = width * (1:nBases) - width/2;
-    for k = 1:nBases
-        idx = ttb(:, k) > ceil(width * (k-1)) & ttb(:, k) <= ceil(width * k);
-        BBstm(idx, k) = 1 / sum(idx);
-    end
-else
-    error('Unknown basis shape');
+switch shape
+    case 'raised cosine'
+        %   ^
+        %  / \
+        % /   \______
+        %      ^
+        %     / \
+        % ___/   \___
+        %         ^
+        %        / \
+        % ______/   \
+        % For raised cosine, the spacing between the centers must be 1/4 of the
+        % width of the cosine
+        dbcenter = nkbins / (3 + nBases); % spacing between bumps
+        width = 4 * dbcenter; % width of each bump
+        bcenters = 2 * dbcenter + dbcenter*(0:nBases-1);
+        % location of each bump centers
+        bfun = @(x,period)((abs(x/period)<0.5).*(cos(x*2*pi/period)*.5+.5));
+        BBstm = bfun(ttb-repmat(bcenters,nkbins,1), width);
+    case 'boxcar'
+        width = nkbins / nBases;
+        BBstm = zeros(size(ttb));
+        bcenters = width * (1:nBases) - width/2;
+        for k = 1:nBases
+            idx = ttb(:, k) > ceil(width * (k-1)) & ttb(:, k) <= ceil(width * k);
+            BBstm(idx, k) = 1 / sum(idx);
+        end
+    otherwise
+        error('Unknown basis shape');
+end
+
+if p.Results.Normalize
+    BBstm = bsxfun(@rdivide, BBstm, sum(BBstm));
+end
+
+if p.Results.Orthogonalize
+    BBstm = orth(BBstm);
 end
 
 bases.type = [shape '@' mfilename];

+buildGLM/addBiasColumn.m

@@ -1,10 +1,14 @@
-function dm = addBiasColumn(dm)
-% Add a column of ones as the first column to estimate the bias (DC term)
+function dm = addBiasColumn(dm, flag)
+% Add a column of ones as the first (or last) column to estimate the bias 
+% (DC term)
 % dm = addBiasColumn(dm);
 %
 % Some regression packages do not allow separate bias estimation, and
 % a constant column of ones must be added to the design matrix. The weight
 % that corresponds to this column will be the bias later.
+%
+% To put the bias column on the right of the design matrix, call:
+% dm = addBiasColumn(dm, 'right')
 
 if nargout ~= 1
     error('Must assign output back to a design matrix!');
@@ -15,5 +19,17 @@
     return
 end
 
-dm.X = [ones(size(dm.X, 1), 1), dm.X];
-dm.biasCol = 1; % indicating that the first column is the bias column
+if nargin < 2
+    flag = 'left'; % bias column defaults to the left of the design matrix
+end
+
+switch flag
+    case 'left'
+        dm.X = [ones(size(dm.X, 1), 1), dm.X];
+        dm.biasCol = 1; % indicating that the first column is the bias column
+    case 'right'
+        n = dm.dspec.edim+1;
+        dm.X = [dm.X, ones(size(dm.X, 1), 1)];
+        dm.biasCol = n; % indicating that the last column is the bias column
+end
+

+buildGLM/addCovariate.m

@@ -64,7 +64,7 @@
 end
 
 if nargin >= 6
-    dspec.covar(newIdx).offset = offset;
+    dspec.covar(newIdx).offset = dspec.expt.binfun(offset);
 else
     dspec.covar(newIdx).offset = 0;
 end

+buildGLM/addCovariateSpiketrain.m

@@ -33,8 +33,7 @@
 assert(ischar(desc), 'Description must be a string');
 
 offset = 1; % Make sure to be causal. No instantaneous interaction allowed.
-
 binfun = dspec.expt.binfun;
-stimHandle = @(trial, expt) basisFactory.deltaStim(binfun(trial.(stimLabel)), binfun(trial.duration));
+stimHandle = @(trial, expt) basisFactory.deltaStim(binfun(trial.(stimLabel)+expt.binSize), binfun(trial.duration));
 
 dspec = buildGLM.addCovariate(dspec, covLabel, desc, stimHandle, basisStruct, offset, varargin{:});

+buildGLM/combineWeights.m

@@ -11,12 +11,20 @@
 
 dspec = dm.dspec;
 binSize = dspec.expt.binSize;
+wout = struct();
 
 if isfield(dm, 'biasCol') % undo z-score operation
-    wout.bias = w(dm.biasCol);
+    if isfield(dm, 'zscore') && numel(dm.zscore.mu) == numel(w) % remove bias from zscore
+        zmu  = dm.zscore.sigma(dm.biasCol);
+        zsig = dm.zscore.mu(dm.biasCol);
+        dm.zscore.sigma(dm.biasCol) = [];
+        dm.zscore.mu(dm.biasCol) = [];
+    else
+        zmu  = 0;
+        zsig = 1;
+    end
+    wout.bias = w(dm.biasCol)*zsig + zmu; % un-z-transform the bias
     w(dm.biasCol) = [];
-    dm.zscore.sigma(dm.biasCol) = [];
-    dm.zscore.mu(dm.biasCol) = [];
 end
 
 if isfield(dm, 'zscore') % undo z-score operation
@@ -35,7 +43,6 @@
 end
 
 startIdx = [1 (cumsum([dspec.covar(:).edim]) + 1)];
-wout = struct();
 
 for kCov = 1:numel(dspec.covar)
     covar = dspec.covar(kCov);

+buildGLM/compileDenseDesignMatrix.m

@@ -8,10 +8,9 @@
 totalT = sum(trialT);
 X      = zeros(totalT, dspec.edim);
 
-trialIndices = trialIndices(:)';
-
-for kTrial = trialIndices
-    ndx = sum(trialT(1:kTrial))-(trialT(kTrial)-1):sum(trialT(1:kTrial));
+for k = 1:numel(trialIndices)
+	kTrial = trialIndices(k);
+    ndx = sum(trialT(1:k))-(trialT(k)-1):sum(trialT(1:k));
 
     for kCov = 1:numel(dspec.covar) % for each covariate
         covar = dspec.covar(kCov);

+buildGLM/getBinnedSpikeTrain.m

@@ -6,11 +6,12 @@
 endTrialIndices = [0 cumsum(binfun([expt.trial(trialIdx).duration])) + 1];
 nT = endTrialIndices(end) - 1; % how many bins total?
 
-for kTrial = trialIdx(:)'
-    bst = endTrialIndices(kTrial) + binfun(expt.trial(kTrial).(spLabel));
-    sts{kTrial} = bst(:);
+for k = 1:numel(trialIdx)
+    kTrial = trialIdx(k);
+    bst = endTrialIndices(k) + binfun(expt.trial(kTrial).(spLabel));
+    sts{k} = bst(:);
 end
 
 sts = cell2mat(sts);
-
+sts(sts>nT) = [];
 y = sparse(sts, 1, 1, nT, 1);