/*
 *  test_random_binomial.sli
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

/* BeginDocumentation
Name: testsuite::test_random_binomial - test of binomial random deviates 

Synopsis: (test_random_binomial) run -> compare binomial random numbers with desired distribution.

Description:
 
  Kolmogorov-Smirnov test of binomial random deviate generator.
  This test performs a Kolmogorov-Smirnov test of the random numbers generated
  by the binomial random deviate generator. It draws realizations and compares 
  the observed cumulative distribution function against the theoretical one.

Remarks:
  The test was implemented following
  http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
  

Author:  July 2011, Moritz Deger
SeeAlso: rdevdict::binomial
*/


/unittest (6688) require
/unittest using

% The factorial is needed to compute the binomial coefficients that occur
% in the probability mass function.
% Direct computation of factorial. If n exceeds 12, integer values can
% not represent the number. This makes the test fail.
/factorial{
  dup 12 leq assert_or_die
  /x Set 1 /ret Set 
  {x 1 lt {exit} if 
  x ret mul /ret Set 
  x 1 sub /x Set } loop 
  ret 
} def


% binomial prob of (k, p, n), k out of n with prob p 
/binomialprob{ /n Set  /p Set  /k Set
  n factorial k factorial  n k sub factorial mul cvd div    %binomial coefficient
  p k pow mul
  1 p sub  n k sub pow mul
} def


% binomial CDF of (k, p, n), k out of n with prob p 
/binomialcdf{
  /n Set /p Set /kl Set
  0.0 /res Set
  0 /cou Set
  {cou kl gt {exit} if           % check if counter > k, exit if true
   cou p n binomialprob res add /res Set    % add probability of stack to return value
   cou 1 add /cou Set                             % increment the counter
   } loop
  res                       % return res
} def


% estimate CDF of observed data
/observedcdf{
  /data Set
  1 add /kp1 Set
  0 /cdfcount Set
  data { kp1 lt {cdfcount 1 add /cdfcount Set} if } forall 
  cdfcount data length_a cvd div
} def


% compute kalpha to compare against 
% python code to get value:
% import scipy.stats
% import scipy.optimize
% alpha = 0.05
% kalpha = scipy.optimize.fmin( lambda x: abs(1-alpha-scipy.stats.ksprob(x)), 1)
0.51962891 /kalpha Set


% actual test, seed as argument
% n should not exceed 12 because then factorial fails
/testbino{
  % get seeded random deviate generator
  /seedn Set
  /n Set
  /p Set
  10000 /nsample Set 
  rngdict /MT19937 get seedn CreateRNG /rng Set
  rng rdevdict /binomial get CreateRDV /bino Set
  
  % set parameters and draw random numbers
  bino << /p p /n n >> SetStatus
  bino nsample RandomArray /data Set

  % test the realizations against the theorical distribution
  0.0 /teststat Set
  0
  { 
  dup n gt {exit} if          % check if stack > n, exit if true 
  dup dup data observedcdf    % 2x counter and observed CDF value on stack
  exch p n binomialcdf        % counter, obs. and theoretical CDF on stack
  
  % compute difference of obs. and theo., memorize only if larger than old value.
  sub dup abs teststat geq {/teststat Set} {pop} ifelse
  1 add                       % increment stack counter
  } loop
  
  % test result (supremum of errors against kolmogorov distribution)
  n sqrt teststat mul kalpha leq assert_or_die
} def


% run the test with several seeds, n and p
0.2 12 11652 testbino
0.1 8 7562 testbino
0.7 11 8766576 testbino