function [k nlin xb STA pSTE pSE] = LNPm_sta_estimate(STE,SE,nB)

% [k nlin xb STA pSTE pSE] = LNPm_sta_estimate(STE,SE,nB)
%
% STE - spike-triggered ensembles: dimensions x ensembes
% SE  - stimulus ensembles: dimensions x ensembes
% nB  - number of bins for non-linearity estimation
%
% k - filter kernel
% nlin - non-parametric non-linearity
% xb - linear response output related to nlin (i.e., the x-axis)
% STA - spike-triggered average
% pSTE, pSE - probability functions for p(s|spike) and p(s)
% 
% Reference:
% Schwartz et al. (2006) Spike-triggered neural characterization. J Vison 6:484-507.
%
% Description: The script estimates a linear filter based on the spike-triggered 
% average (STA) and estimates a non-parametric non-linearity. Hence, the script 
% provides linear-nonlinear model estimates.
%
% ---------
%
%    Copyright (C) 2016, B. Herrmann
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% ----------------------------------------------------------------------------
% B. Herrmann, Email: herrmann.b@gmail.com, 2016-04-11

% spike-triggered average
STA = mean(STE,2);

% filter kernel
k = STA - mean(STA);

% dot product / STA response / linear filter output
dprodSTE = k' * STE;
dprodSE  = k' * SE;

% bin linear output
rB  = [-max(abs(dprodSE)) max(abs(dprodSE))];
bSE = linBin(dprodSE,rB,diff(rB)/nB,nB);
[bSTE xb] = linBin(dprodSTE,rB,diff(rB)/nB,nB);

% get proportion/probability function
pSE  = bSE/length(dprodSE);
pSTE = bSTE/length(dprodSTE);

% non-parametric estimation of the non-linearity
nlin = pSTE./pSE;



% -------------------------
% sub functions
% -------------------------

function [B xb] = linBin(D,rB,width,nB)
% D - data vector; rB - overall bounds; width - width of bins; nB - number of bins
% B - binned data; xb - x-axis

% get bins
xb    = linspace(rB(1)+width/2,rB(2)-width/2,nB);   % center of bins: lin
xbins = [xb'-width/2 xb'+width/2];                  % edges of bins: lin

B = zeros([size(xbins,1) 1]);
for ii = 1 : size(xbins,1)
	idxs = D >= xbins(ii,1) & D <= xbins(ii,2);
	B(ii) = length(D(idxs));
end

