#! /usr/bin/env python
#
def ignnbn ( n, p ):

#*****************************************************************************80
#
## IGNNBN generates a negative binomial random deviate.
#
#  Discussion:
#
#    This procedure generates a single random deviate from a negative binomial
#    distribution.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    Original FORTRAN77 version by Barry Brown, James Lovato.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Luc Devroye,
#    Non-Uniform Random Variate Generation,
#    Springer, 1986,
#    ISBN: 0387963057,
#    LC: QA274.D48.
#
#  Parameters:
#
#    Input, integer N, the required number of events.
#    0 <= N.
#
#    Input, real P, the probability of an event during a
#    Bernoulli trial.  0.0 < P < 1.0.
#
#    Output, integer VALUE, a random deviate from
#    the distribution.
#
  from gengam import gengam
  from ignpoi import ignpoi
  from sys import exit

  if ( n < 0 ):
    print ( '' )
    print ( 'IGNNBN - Fatal error!' )
    print ( '  N < 0.' )
    exit ( 'IGNNBN - Fatal error!' )

  if ( p <= 0.0 ):
    print ( '' )
    print ( 'IGNNBN - Fatal error!' )
    print ( '  P <= 0.0' )
    exit ( 'IGNNBN - Fatal error!' )

  if ( 1.0 <= p ):
    print ( '' )
    print ( 'IGNNBN - Fatal error!' )
    print ( '  1.0 <= P' )
    exit ( 'IGNNBN - Fatal error!' )
#
#  Generate Y, a random gamma (n,(1-p)/p) variable.
#
  r = float ( n )
  a = p / ( 1.0 - p )
  y = gengam ( a, r )
#
#  Generate a random Poisson ( y ) variable.
#
  value = ignpoi ( y )

  return value

def ignnbn_test ( phrase ):

#*****************************************************************************80
#
## IGNNBN_TEST tests IGNNBN, which generates Negative Binomial deviates.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    03 September 2018
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from genunf import genunf
  from initialize import initialize
  from phrtsd import phrtsd
  from set_initial_seed import set_initial_seed
  from stats import stats
  from trstat import trstat

  n = 10000

  print ( '' )
  print ( 'IGNNBN_TEST' )
  print ( '  IGNNBN generates negative binomial deviates.' )
#
#  Initialize the generators.
#
  initialize ( )
#
#  Set the seeds based on the phrase.
#
  seed1, seed2 = phrtsd ( phrase )
#
#  Initialize all generators.
#
  set_initial_seed ( seed1, seed2 )
#
#  Select the parameters at random within a given range.
#
  low = 3.0
  high = 20.0
  nn = int ( genunf ( low, high ) )

  low = 0.0
  high = 1.0
  pp = genunf ( low, high )

  print ( '' )
  print ( '  N = %d' % ( n ) )
  print ( '' )
  print ( '  Parameters:' )
  print ( '' )
  print ( '  NN = %d' % ( nn ) )
  print ( '  PP = %g' % ( pp ) )
#
#  Generate N samples.
#
  array = np.zeros ( n )
  for i in range ( 0, n ):
    array[i] = ignnbn ( nn, pp )
#
#  Compute statistics on the samples.
#
  av, var, xmin, xmax = stats ( array, n )
#
#  Request expected value of statistics for this distribution.
#
  pdf = 'nbn'
  param = np.array ( [ nn, pp ] )
 
  avtr, vartr = trstat ( pdf, param )

  print ( '' )
  print ( '  Sample data range:          %14.6g  %14.6g' % ( xmin, xmax ) )
  print ( '  Sample mean, variance:      %14.6g  %14.6g' % ( av,   var ) )
  print ( '  Distribution mean, variance %14.6g  %14.6g' % ( avtr, vartr ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'IGNNBN_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  phrase = 'Randomizer'
  ignnbn_test ( phrase )
  timestamp ( )