function PDF = BinomPDF(model,v);
% PDF = PDF(model,v);
% Probability density function of simple binomial with noise, quantal variance and stim failurenv=length(v);
% NB v has to be an array: [ 1 nv ] (and not [nv 1]).

if model.p<0 || model.p>1 || model.pstim<0 || model.pstim>1 || model.q<0 || model.qsigma<0 || model.noise<0 ...
        || abs(model.offset) > 0.5*model.q 
   PDF = 0;
else
   % Remove offset:
   v = v - model.offset;
   
   nv=length(v);
   PDF=zeros(1,nv);
   probterms=zeros(1,model.n+1);
   m=[0:model.n];
   for j=1:model.n+1
%      probterms(j) =  prod([model.n-m(j)+1:model.n]) / prod([1:m(j)]) * model.p^m(j) * (1-model.p)^(model.n-m(j)) ;
      probterms(j) =  nchoosek(model.n,m(j)) * model.p^m(j) * (1-model.p)^(model.n-m(j)) ;
   end    
   v2d = repmat(v,[model.n+1 1]);
   m2d = repmat(m',[1 nv]);
   probterms2d = repmat(probterms',[1 nv]);
   
   if (model.QVtype==1)
      sd = sqrt(model.noise.^2+m.*model.qsigma.^2);
      sd2d = repmat(sd',[1 nv]);      
   elseif model.QVtype==2 % "flat"
      % SD = noise for zeroth peak, = (noise+qsig) for subsequent peaks
      var = model.noise.^2+(m>0).*model.qsigma.^2;
      sd = sqrt(var);
      sd2d = repmat(sd',[1 nv]);
   end
   PDF =  model.pstim * sum(probterms2d .* exp( -(v2d-model.q.*m2d).^2./2./sd2d.^2 ) ./ sd2d )...
      + (1-model.pstim) * exp( -v.^2./2./model.noise.^2 ) / model.noise;
   PDF = PDF ./ sqrt(2*pi);
   
   % 15th May 2008
   % Now add in "glitch rate"
   % There is a finite chance that instead of being generated by a quantal
   % event + noise as above, an EPSP is just a fluke. Assume that glitches
   % are evenly distributed between model.glitchmin and model.glitchmax
   if isfield(model,'glitchrate')
       PDF = (1-model.glitchrate) * PDF + model.glitchrate / (model.glitchmax-model.glitchmin) * ( v>model.glitchmin & v<model.glitchmax );
   end   
end
