getSurprise

  getSurprise method for jpeth class

  Examples:
  surprise=getSurprise(obj)
  [surprise excite inhib]=getSurprise(obj)

  Returns the "surprise" matrix for the jpeth object as defined in Palm et
  al. 1988:
            surprise=excite-inhib
  where
            excite=-ln(prob(z>=M given K and L)) and
            inhib=-ln(prob(z<=M given K and L)
  M is number of coincidences in the jpeth raw matrix and K and L are the
  corresponding number of spikes in peth1 and peth2 of the jpeth object
  (reversed if K>L). Corrections are included for Palm's Type 1 and 2
  errors following his original FORTRAN code.

  Algorithm:
  The probabilities of z==M are estimated following Eqn. 7 of Palm et
  al. 1988 but natural logs are used throughout to minimize problems due
  to IEEE precision and to enhance speed: 
            ln(x!) is estimated as gammaln(x+1) 
            ln(choosek(n,k)) as gammaln(n+1)-(gammaln(n-k+1)+gammaln(k+1));
                (see Press et al. 1992)
  Eqn 7 is then evaluated by log addition/subtraction before finally taking
  the exponent to return probabilities. This avoids intermediate results
  beyond the 15 significant digits of IEEE 64-bit floating point
  arithmetic

  References
  Palm et al. (1988) On the significance of correlations among neural spike
  trains Biological Cybernetics 59, 1-11.
  Press et al. (1992) Numerical Recipes in C.


  See also jpeth

  -------------------------------------------------------------------------
  Author: Malcolm Lidierth 01/09
  Copyright © The Author & King's College London 2009-
  -------------------------------------------------------------------------