[CONTACT]

[ABOUT]

[POLICY]

For other Languages simply replace

Found at: ftp.icm.edu.pl:70/packages/netlib/random/amrandom.pas

{ Resource Strings for use with AMRandom.

	For other Languages simply replace this unit.}

unit AMRandomRS;

nterface

ResourceString

	rsDivideByZero = 'Divide By Zero';
	rsInvalidShape = 'Invalid Shape';
	rsShapeTooSmall = 'Shape too small';
	rsInvalidValue = 'Invalid Value';
	rsVectorEmpty = 'Vector Empty';
	rsVectorTooSmall = 'Vector too small';

mplementation

end.




{ Delphi Translation developed by ESB Consultancy
	http://www.esbconsult.com.au/
	mailto:info@esbconsult.com.au

	v1.1.2 9-Jan-2002
	- Borland C++ Builder Support added

	v1.1.1 19-Sep-2001
	- Fixed bug in Multivariate Normal Distribution - thanks to Haris Baltzakis

	v1.1 19-Sep-2001
	- now tested with Delphi 6
	- Tester now built with Delphi 6 and ESBPCS v2

	v1.00 9-Jan-2001
	- Fixed a bug that Alan Miller reported in Random_T
	- Some Improvements to the TestRandom Project
	- Now also incorporated into our ESBPCS with the release of v1.5

	Public Beta 3 Release 23-Oct-2000
	- Fixed problem where Poisson would get divide by zero - incorrect initialisation
		on my part.

	Public Beta 2 Release 12-Oct-2000
	- Exit condition for Random_Normal modified slightly
	- Changed routines to be overloaded to allow you to pass in the Uniform
		Random Number Generator you want to use. If omitted then Delphi's
		Random Number Generator will be used.
	- Added routines to call the Delphi Random and the MRNG Random Generators
		in the form of TRandomGenFunction.
	- Added AMRandomRS Unit to store strings as Resource Strings. Replace this
		unit with other language translations. Feel free to send to us so
		we can share with all users.
	- Added a Test Utility - source is included but ESBPCS (Shareware) and
		EFD Systems' MRNG.PAS (Freeware) is required.
		ESBPCS: http://www.esbconsult.com.au/esbpcs.html
		EFD: http://efd.home.mindspring.com/tools.htm
	- Test Utility still doesn't test Multi-variate.

	Public Beta 1: 30-Sep-2000
}

{ A unit for random number generation from the following distributions:

	Distribution                    Function/subroutine name

	Normal (Gaussian)               Random_Normal
	Gamma                           Random_Gamma
	Chi-squared                     Random_ChiSq
	Exponential                     Random_Exponential
	Weibull                         Random_Weibull
	Beta                            Random_Beta
	t                               Random_T
	Multivariate normal             Random_MVNorm
	Generalized inverse Gaussian    Random_Inv_Gauss
	Poisson                         Random_Poisson
	Binomial                        Random_Binomial1   *
							  Random_Binomial2   *
	Negative binomial               Random_neg_Binomial
	von Mises                       Random_von_Mises
	Cauchy                          Random_Cauchy

	Generate a random ordering of the integers 1 .. N
							  Random_Order
	Initialize (seed) the uniform random number generator for ANY compiler
							  Seed_Rrandom_Number

  * Two functions are provided for the binomial distribution.
  If the parameter values remain constant, it is recommended that the
  first function is used (random_binomial1).   If one or both of the
  parameters change, use the second function (random_binomial2).

  Delphi's own random number generator, Random, is used to provide a source of
  uniformly distributed random numbers.

  N.B. At this stage, only one random number is generated at each call to
	  one of the functions above.

 The unit uses the following functions which are included here:
 bin_prob to calculate a single binomial probability
 lngamma  to calculate the logarithm to base e of the gamma function

 Some of the code is adapted from Dagpunar's book:
	Dagpunar, J. 'Principles of random variate generation'
	Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 In most of Dagpunar's routines, there is a test to see whether the value
 of one or two floating-point parameters has changed since the last call.
 These tests have been replaced by using a logical variable FIRST.
 This should be set to .TRUE. on the first call using new values of the
 parameters, and .FALSE. if the parameter values are the same as for the
 previous call.

 First Delphi Version by Glenn Crouch of ESB Consultancy
					info@esbconsult.com.au
					http://www.esbconsult.com.au/

 FORTRAN Version History
 Version 1.12, 2 July 2000
 Changes from version 1.01
 1. The random_order, random_Poisson & random_binomial routines have been
    replaced with more efficient routines.
 2. A routine, seed_random_number, has been added to seed the uniform random
    number generator.   This requires input of the required number of seeds
    for the particular compiler from a specified I/O unit such as a keyboard.
 3. Made compatible with Lahey's ELF90.
 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.

	Author: Alan Miller
		   CSIRO Division of Mathematical & Information Sciences
		   Private Bag 10, Clayton South MDC
		   Clayton 3169, Victoria, Australia
	Phone: (+61) 3 9545-8036      Fax: (+61) 3 9545-8080
	e-mail: Alan.Miller @ vic.cmis.csiro.au

}
unit AMRandom;

nterface

{.$Define UseMRNG} // Remove this if you don't have EFD's free MRNG.PAS

type
	TAMFloatVector = array of Extended;

type
	TRandomGenFunction = function: Extended;

{: Used to call Delphi's Inbuilt Random Number Generator. Remember to use
	Randomize if you don't want repeatable sequences. }
function DelphiRandom: Extended;

{$IFDEF UseMRNG}
{: Used to call EFD System's MRNG Random Number Generator. }
function MRNGRandom: Extended;
{$ENDIF}

{ Adapted from the following Fortran 77 code
	 ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
	 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
	 VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

  The function random_normal() returns a normally distributed pseudo-random
  number with zero mean and unit variance.

  The algorithm uses the ratio of uniforms method of A.J. Kinderman
  and J.F. Monahan augmented with quadratic bounding curves. }
function Random_Normal: Extended; overload;
function Random_Normal (RandomGenerator: TRandomGenFunction): Extended; overload;


{ Adapted from Fortran 77 code from the book:
	Dagpunar, J. 'Principles of random variate generation'
	Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

	FUNCTION GENERATES A RANDOM GAMMA VARIATE.
	CALLS EITHER random_gamma1 (Shape > 1.0)
	OR random_exponential (Shape = 1.0)
	OR random_gamma2 (Shape < 1.0).

	Shape = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
}
function Random_Gamma (const Shape: Extended): Extended; overload;
function Random_Gamma (const Shape: Extended;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{  	Generates a random variate from the chi-squared distribution with
	dff degrees of freedom.
}
function Random_ChiSq (const DF: Integer): Extended; overload;
function Random_ChiSq (const DF: Integer;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
  TO EXP(-random_exponential), USING INVERSION.
}
function Random_Exponential: Extended; overload;
function Random_Exponential (RandomGenerator: TRandomGenFunction): Extended; overload;

{     Generates a random variate from the Weibull distribution with
	probability density:
				  a
			a-1  -x
	f(x) = a.x    e
}
function Random_Weibull (const a: Extended): Extended; overload;
function Random_Weibull (const a: Extended;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
  FROM A BETA DISTRIBUTION WITH DENSITY
  PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
  USING CHENG'S LOG LOGISTIC METHOD.

	 AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
	 BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
}

function Random_Beta (const aa, bb: Extended): Extended; overload;
function Random_Beta (const aa, bb: Extended;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{  Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE FROM A
  T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.

	 DF = DEGREES OF FREEDOM OF DISTRIBUTION
		  (DF >= 1)
}
function Random_T (const DF: Integer): Extended; overload;
function Random_T (const DF: Integer;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{  Adapted from Fortran 77 code from the book:
	Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  N.B. An extra argument, ier, has been added to Dagpunar's routine

	 SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
	 VECTOR USING A CHOLESKY DECOMPOSITION.

  ARGUMENTS:
	    N = NUMBER OF VARIATES IN VECTOR
		  (INPUT,INTEGER >= 1)
	 H(J) = J'TH ELEMENT OF VECTOR OF MEANS
		  (INPUT,REAL)
	 X(J) = J'TH ELEMENT OF DELIVERED VECTOR
		  (OUTPUT,REAL)

	D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
		   (INPUT,REAL)
	F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
		 DECOMPOSITION OF VARIANCE MATRIX (J <= I)
		   (OUTPUT,REAL)

    ier = 1 if the input covariance matrix is not +ve definite
	   = 0 otherwise
}
	var f, x: TAMFloatVector; var ier: Integer); overload;
	var f, x: TAMFloatVector; var ier: Integer;
	RandomGenerator: TRandomGenFunction); overload;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM
  A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
  WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
  USING A RATIO METHOD.

	 h = PARAMETER OF DISTRIBUTION (0 <= REAL)
	 b = PARAMETER OF DISTRIBUTION (0 < REAL)
}

function Random_Inv_Gauss (const h, b: Extended): Extended; overload;
function Random_Inv_Gauss (const h, b: Extended;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{  FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
  This algorithm is suitable when many random variates are required
  with the SAME parameter values for n & p.

	P = BERNOULLI SUCCESS PROBABILITY
		  (0 <= REAL <= 1)
	N = NUMBER OF BERNOULLI TRIALS
		  (1 <= INTEGER)

  Reference: Kemp, C.D. (1986). `A modal method for generating binomial
		   variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
}
function  Random_Binomial1 (const n: Integer; const p: Extended): Int64; overload;
function  Random_Binomial1 (const n: Integer; const p: Extended;
	RandomGenerator: TRandomGenFunction): Int64; overload;

{**********************************************************************
	 Translated to Fortran 90 by Alan Miller from:
						 RANLIB

	 Library of Fortran Routines for Random Number Generation

				   Compiled and Written by:

					   Barry W. Brown
					    James Lovato

			 Department of Biomathematics, Box 237
			 The University of Texas, M.D. Anderson Cancer Center
			 1515 Holcombe Boulevard
			 Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

				 GENerate BINomial random deviate

						 Function

	 Generates a single random deviate from a binomial
	 distribution whose number of trials is N and whose
	 probability of an event in each trial is P.

						 Arguments

	 N  --> The number of trials in the binomial distribution
		   from which a random deviate is to be generated.
						 INTEGER N

	 P  --> The probability of an event in each trial of the
		   binomial distribution from which a random deviate
		   is to be generated.
						 REAL P

}
function Random_Binomial2 (const n: Int64; const pp: Extended): Int64; overload;
function Random_Binomial2 (const n: Int64; const pp: Extended;
	RandomGenerator: TRandomGenFunction): Int64; overload;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
  INVERSION AND/OR THE REPRODUCTIVE PROPERTY.

	SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
	   = the `power' parameter of the negative binomial
		  (0 < REAL)
	P = BERNOULLI SUCCESS PROBABILITY
		  (0 < REAL < 1)

  THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
  OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
  THE REPRODUCTIVE PROPERTY IS USED.
}
function Random_Neg_Binomial (const sk, p: Extended): Int64; overload;
function Random_Neg_Binomial (const sk, p: Extended;
	RandomGenerator: TRandomGenFunction): Int64; overload;

{     Algorithm VMD from:
	 Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
	 comparison of random numbers', J. of Appl. Statist., 17, 165-168.

	 Fortran 90 code by Alan Miller
	 CSIRO Division of Mathematical & Information Sciences

	 Arguments:
	 k (real)        parameter of the von Mises distribution.
}

function Random_von_Mises (const k: Extended): Extended; overload;
function Random_von_Mises (const k: Extended;
	RandomGenerator: TRandomGenFunction): Extended; overload;

{      Generate a random deviate from the standard Cauchy distribution
}
function Random_Cauchy: Extended; overload;
function Random_Cauchy (RandomGenerator: TRandomGenFunction): Extended; overload;

{**********************************************************************
	 Translated to Fortran 90 by Alan Miller from:
					   RANLIB

	 Library of Fortran Routines for Random Number Generation

				 Compiled and Written by:

					 Barry W. Brown
					  James Lovato

		    Department of Biomathematics, Box 237
		    The University of Texas, M.D. Anderson Cancer Center
		    1515 Holcombe Boulevard
		    Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

				 GENerate POIsson random deviate

					    Function

  Generates a single random deviate from a Poisson distribution with mean mu.

					    Arguments

	 mu --> The mean of the Poisson distribution from which
		   a random deviate is to be generated.
						 REAL mu

						 Method

	 For details see:

			 Ahrens, J.H. and Dieter, U.
			 Computer Generation of Poisson Deviates
			 From Modified Normal Distributions.
			 ACM Trans. Math. Software, 8, 2
			 (June 1982),163-179

	 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
	 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
}
function Random_Poisson(var mu: Extended): Int64; overload;
function Random_Poisson(var mu: Extended;
	RandomGenerator: TRandomGenFunction): Int64; overload;

{ Logarithm to base e of the gamma function.

  Accurate to about 1.e-14.
  Programmer: Alan Miller

  Latest revision of Fortran 77 version - 28 February 1988
}
function lngamma (const x: Extended): Extended;

mplementation

uses
	{$IFDEF UseMRNG}
	MRNG,
	{$ENDIF}
	Math, SysUtils,
	AMRandomRS;

const
	VSmall = MinDouble;
	VLarge = MaxDouble;

function DelphiRandom: Extended;
begin
	Result := Random;
end;

{$IFDEF UseMRNG}
function MRNGRandom: Extended;
begin
	Result := MRandom;
end;
{$ENDIF}


{ Adapted from the following Fortran 77 code
	 ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
	 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
	 VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

  The function random_normal() returns a normally distributed pseudo-random
  number with zero mean and unit variance.

  The algorithm uses the ratio of uniforms method of A.J. Kinderman
  and J.F. Monahan augmented with quadratic bounding curves. }
function Random_Normal (RandomGenerator: TRandomGenFunction): Extended;
const
	s = 0.449871;
	t = -0.386595;
	a = 0.19600;
	b = 0.25472;
	r1 = 0.27597;
	r2 = 0.27846;
var
	u, v, x, y, q: Extended;
	Done: Boolean;

begin
//     Generate P = (u,v) uniform in rectangle enclosing acceptance region

	Done := False;
	repeat
		u := RandomGenerator;
		v := RandomGenerator;
		v := 1.7156 * (v - 0.5);

		// Evaluate the quadratic form
		x := u - s;
		Y := abs (v) - t;
		q := Sqr (x) + y * (a * y - b * x);

		// Accept P if inside inner ellipse
		if (q < r1) then
			Done := True
		else if (q <= r2) and (Sqr (v) < -4.0 * Ln (u) * Sqr (u)) then
			Done := True;
	until Done;
	// Return ratio of P's coordinates as the normal deviate
	if u < VSmall then
		raise EMathError.Create (rsDivideByZero);

	Result := v / u;
end;

function Random_Normal: Extended;
begin
	Result := Random_Normal (DelphiRandom);
end;

{ Uses the algorithm in
  Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
  gamma variables', Trans. om Math. Software (TOMS), vol.26.

  Generates a random gamma deviate for shape parameter Shape >= 1.
}
function Random_Gamma1 (const Shape: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
var
	c, d: Extended;
	u, v, x: Extended;
	Done: Boolean;
begin
	if (Shape <= 1.0) then
		raise EMathError.Create (rsInvalidShape);

	d := Shape - 1 / 3;
	c := 1 / SQRT (9.0 * d);

	Result := 0.0;

	Done := False;
	repeat
		// Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.

		repeat
			x := Random_Normal (RandomGenerator);
			v := Power (1 + c * x, 3);
		until v > 0;

		// Generate uniform variable U

		u := RandomGenerator;

		if (u < 1.0 - 0.0331 * Sqr (Sqr (x))) then
		begin
			Result := d * v;
			Done := True;
		end
		else if Ln (u) < 0.5 * Sqr (x) + d * (1.0 - v + Ln (v)) then
		begin
			Result := d * v;
			Done := True;
		end
	until Done;
end;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
  GAMMA2**(S-1) * EXP(-GAMMA2),
  USING A SWITCHING METHOD.

	S = SHAPE PARAMETER OF DISTRIBUTION
		 (REAL < 1.0)
 }
function Random_Gamma2 (const Shape: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
var
	r, x, w: Extended;
	a, p, c, uf, vr, d: Extended;
	Done: Boolean;
begin
	if (Shape <= 0.0) and (Shape >= 1.0) THEN
		raise EMathError.Create (rsInvalidShape);
	if (Shape < VSmall) then
		raise EMathError.Create (rsShapeTooSmall);

	a := 1.0 - Shape;
	p := a / (a + Shape * EXP(-a));
	c := 1.0 / Shape;
	uf := p * Power (VSmall / a, Shape);
	vr := 1.0 - VSmall;
	d := a * Ln (a);
	w := 0;
	x := 0;

	Done := False;
	repeat
		r := RandomGenerator;
		if (r >= vr) then
			Continue;

		if (r > p) then
		begin
			x := a - Ln ((1.0 - r)/(1.0 - p));
			w := a * Ln (x) - d;
		end
		else if (r > uf) then
		begin
			x := a * Power (r / p, c);
			w := x;
		end
		else
		begin
		    Result := 0.0;
		    Exit;
		end;

		r := RandomGenerator;

		if (1.0 - r <= w) and (r > 0.0) then
		begin
			if (r * (w + 1.0) >= 1.0 ) then
				Continue;

			if (-Ln (r) <= w)  then
				Continue;

			Done := True;
		end;
	until Done;
	Result := x;
end;

{ Adapted from Fortran 77 code from the book:
	Dagpunar, J. 'Principles of random variate generation'
	Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

	FUNCTION GENERATES A RANDOM GAMMA VARIATE.
	CALLS EITHER random_gamma1 (Shape > 1.0)
	OR random_exponential (Shape = 1.0)
	OR random_gamma2 (Shape < 1.0).

	Shape = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
}
function Random_Gamma (const Shape: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
begin
	if (Shape <= 0.0) then
		raise EMathError.Create (rsInvalidShape);

	if (Shape > 1.0) then
		Result := Random_Gamma1 (Shape, RandomGenerator)
	else if (Shape < 1.0) then
		Result := Random_Gamma2 (Shape, RandomGenerator)
	else
		Result := Random_Exponential (RandomGenerator);
end;

function Random_Gamma (const Shape: Extended): Extended;
begin
	Result := Random_Gamma (Shape, DelphiRandom);
end;


{  	Generates a random variate from the chi-squared distribution with
	dff degrees of freedom.
}
function Random_ChiSq (const DF: Integer;
	RandomGenerator: TRandomGenFunction): Extended;
begin
	Result := 2.0 * Random_Gamma (0.5 * DF, RandomGenerator)
end;

function Random_ChiSq (const DF: Integer): Extended;
begin
	Result := Random_ChiSq (DF, DelphiRandom);
end;


{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
  TO EXP(-random_exponential), USING INVERSION.
}
function  Random_Exponential (RandomGenerator: TRandomGenFunction): Extended;
var
	r: Extended;
	Done: Boolean;
begin
	Done := False;
	repeat
		r := RandomGenerator;
		if (r > vSmall) then
			Done := True;
	until Done;
	Result := - Ln (r)
end;

function  Random_Exponential : Extended;
begin
	Result := Random_Exponential (DelphiRandom);
end;

{     Generates a random variate from the Weibull distribution with
	probability density:
				  a
			a-1  -x
	f(x) = a.x    e
}
function Random_Weibull (const a: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
begin
	if Abs (a) <= VSmall then
		raise EMathError.Create (rsInvalidValue);

	Result := Power (Random_Exponential (RandomGenerator), 1.0 / a);
end;

function Random_Weibull (const a: Extended): Extended;
begin
	Result := Random_Weibull (a, DelphiRandom);
end;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
  FROM A BETA DISTRIBUTION WITH DENSITY
  PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
  USING CHENG'S LOG LOGISTIC METHOD.

	 AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
	 BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
}

function Random_Beta (const aa, bb: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
const
	aln4 = 1.3862943611198906;
var
	a, b, g, r, s, x, y, z: Extended;
	d, f, h, t, c: Extended;
	Swap: Boolean;
	Done: Boolean;
begin
	if (aa <= 0.0) or (bb <= 0.0) then
		raise EMathError.Create (rsInvalidValue);

	a := aa;
	b := bb;
	swap := b > a;
	if swap then
	begin
		g := b;
		b := a;
		a := g;
	end;
	d := a / b;
	f := a + b;
	if (b > 1.0) then
	begin
		h := Sqrt ((2.0 * a *b - f) / (f - 2.0));
		t := 1.0;
	end
	else
	begin
		h := b;
		t := 1.0 / (1.0 + Power (a / (VLarge * b), b));
	end;
	c := a + h;

	Done := False;
	Result := 0.0;
	repeat
		r := RandomGenerator;
		x := RandomGenerator;
		s := Sqr (r) * x;
		if (r < VSmall) or (s <= 0.0) then
			Continue;

		if (r < t) then
		begin
			x := Ln (r / (1.0 - r)) / h;
			y := d * Exp (x);
			z := c * x + f * Ln ((1.0 + d)/(1.0 + y)) - aln4;
			if (s - 1.0 > z) then
			begin
				if (s - s * z > 1.0) then
					Continue;
				if (Ln (s) > z) then
					Continue;
			end;
			Result := y / (1.0 + y);
		end
		else
		begin
			if 4.0 * s > Power (1.0 + 1.0 / d, f) then
				Continue;
			Result := 1.0
		end;
		Done := True
	until Done;

	if Swap then
		Result := 1 - Result;
end;

function Random_Beta (const aa, bb: Extended): Extended;
begin
	Result := Random_Beta (aa, bb, DelphiRandom);
end;

{  Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE FROM A
  T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.

	 DF = DEGREES OF FREEDOM OF DISTRIBUTION
		  (DF >= 1)
}
function Random_T (const DF: Integer;
	RandomGenerator: TRandomGenFunction): Extended;
var
	s, c, a, f, g: Extended;
	r, x, v: Extended;
	Done: Boolean;
begin
	if (DF < 1) then
		raise EMathError.Create (rsInvalidValue);

	s := DF;
	c := -0.25 * (s + 1.0);
	a := 4.0 / Power (1.0 + 1.0 / s, c);
	f := 16.0 / a;
	if (DF > 1) then
	begin
		g := s - 1.0;
		g := Power ((s + 1.0) /g , c) * Sqrt ((s + s) / g);
	end
	else
		g := 1.0;

	Done := False;
	x := 0;
	repeat
		r := RandomGenerator;
		if (r <= 0.0) then
			Continue;
		v := RandomGenerator;
		x := (2.0 * v - 1.0) * g / r;
		v := Sqr (x);
		if (v > 5.0 - a * r) then
		begin
			if (DF >= 3) and (r * (v + 3.0) > f) then
				Continue;
			if r > Power (1.0 + v / s, c) then
				Continue;
		end;
		Done := True;
	until Done;
	Result := x;
end;

function Random_T (const DF: Integer): Extended;
begin
	Result := Random_T (DF, DelphiRandom);
end;

{  Adapted from Fortran 77 code from the book:
	Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  N.B. An extra argument, ier, has been added to Dagpunar's routine

	 SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
	 VECTOR USING A CHOLESKY DECOMPOSITION.

  ARGUMENTS:
	    N = NUMBER OF VARIATES IN VECTOR
		  (INPUT,INTEGER >= 1)
	 H(J) = J'TH ELEMENT OF VECTOR OF MEANS
		  (INPUT,REAL)
	 X(J) = J'TH ELEMENT OF DELIVERED VECTOR
		  (OUTPUT,REAL)

	D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
		   (INPUT,REAL)
	F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
		 DECOMPOSITION OF VARIANCE MATRIX (J <= I)
		   (OUTPUT,REAL)

    ier = 1 if the input covariance matrix is not +ve definite
	   = 0 otherwise
}
	var f, x: TAMFloatVector; var ier: Integer;
	RandomGenerator: TRandomGenFunction);
var
	i, j, m, n, n2, fn: Integer;
	y, v: Extended;
begin
	n := Length (h);
	if (n = 0) then
		raise EMathError.Create (rsVectorEmpty);

	fn := n * (n + 1) div 2;
	if (Length (D) < fn) then
		raise EMathError.Create (rsVectorTooSmall);

	SetLength (f, fn);

	ier := 0;

	n2 := 2 * n;
	if (d [0] <= 0.0) then
	begin
		ier := 1;
		Exit;
	end;

	f [0] := Sqrt (d [0]);
	y := 1.0 / f [0];

	for j := 2 to n do
		f [j - 1] := d [j * (j - 1) div 2] * y;

	for i := 2 to n do
	begin
		v := d [i * (i - 1) div 2 + i - 1];
		for m := 1 to i - 1 do
			v := v - Sqr (f [(m - 1) *(n2 - m) div 2 +  i - 1]);

		if (v <= 0.0) then
		begin
			ier := 1;
			Exit;
		end;

		v := Sqrt (v);
		y := 1.0 / v;
		F [(i - 1) * (n2 - i) div 2 + i - 1] := v;

		for j := i + 1 to n do
		begin
			v := d [j * (j - 1) div 2 + i - 1];
			for m := 1 to i-1 do
			begin
				v := v - f [(m - 1) *(n2 - m) div 2 + i - 1]
					* f [(m - 1) * (n2 - m) div 2 + j - 1]
			end;
			f [(i - 1) * (n2 - i) div 2 + j - 1] := v * y;
		end;
	end;

	X := Copy (H, 0, n);
	for j := 1 to n do
	begin
		y := Random_Normal (RandomGenerator);
		for i := j to n do
			x [i - 1] := x [i - 1] + f [(j - 1) *(n2 - j) div 2 + i - 1] * y;
	end;
end;

	var f, x: TAMFloatVector; var ier: Integer);
begin
	Random_MVNorm (h, d, f, x, ier, DelphiRandom);
end;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM
  A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
  WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
  USING A RATIO METHOD.

	 h = PARAMETER OF DISTRIBUTION (0 <= REAL)
	 b = PARAMETER OF DISTRIBUTION (0 < REAL)
}

function Random_Inv_Gauss (const h, b: Extended;
	RandomGenerator: TRandomGenFunction): Extended; 
var
	ym, xm, r, w, r1, r2, x: Extended;
	a, c, d, e: Extended;
	Done: Boolean;
begin
	if (h < 0.0) or (b <= 0.0) then
		raise EMathError.Create (rsInvalidValue);

	if (h > 0.25 * b * Sqrt (VLarge)) then
		raise EMathError.Create (rsInvalidValue);

	e := Sqr (b);
	d := h + 1.0;
	ym := (-d + Sqrt (Sqr (d) + e)) / b;
	if (ym < VSmall) then
		raise EMathError.Create (rsInvalidValue);

	d := h - 1.0;
	xm := (d + Sqrt (Sqr (d) + e)) / b;
	if (xm < VSmall) then
		raise EMathError.Create (rsInvalidValue);

	d := 0.5 * d;
	e := -0.25 * b;
	r := xm + 1.0 / xm;
	w := xm * ym;
	a := Power (w, -0.5 *h) * Sqrt (xm / ym) * Exp (-e *(r - ym - 1.0 / ym));
	if (a < vsmall) then
		raise EMathError.Create (rsInvalidValue);

	c := -d * Ln (xm) - e * r;

	Done := False;
	x := 0.0;
	repeat
		r1 := RandomGenerator;
		if (r1 > 0.0) then
		begin
			r2 := RandomGenerator;
			x := a * r2 / r1;
			if (x > 0.0) then
			begin
				if (Ln (r1) < d * Ln (x) + e * (x + 1.0 / x) + c) then
					Done := True
			end;
		end
	until Done;
	Result := x;
end;

function Random_Inv_Gauss (const h, b: Extended): Extended; 
begin
	Result := Random_Inv_Gauss (h, b, DelphiRandom);
end;

{ Logarithm to base e of the gamma function.

  Accurate to about 1.e-14.
  Programmer: Alan Miller

  Latest revision of Fortran 77 version - 28 February 1988
}
function lngamma (const x: Extended): Extended;
const
	a1 = -4.166666666554424E-02;
	a2 = 2.430554511376954E-03;
	a3 = -7.685928044064347E-04;
	a4 = 5.660478426014386E-04;
	lnrt2pi = 9.189385332046727E-1;
	pi = 3.141592653589793E0;
var
	temp, arg, product: Extended;
	reflect: Boolean;
begin
//  lngamma is not defined if x = 0 or a negative integer.
	if (x = 0.0) or ((x < 0.0) and (Abs (x - Int (x)) < VSmall)) then
		raise EMathError.Create (rsInvalidValue);

// If x < 0, use the reflection formula:
//        gamma(x) * gamma(1-x) = pi * cosec(pi.x)

	reflect := x < 0.0;
	if reflect then
		arg := 1.0 - x
	else
		arg := x;

// Increase the argument, if necessary, to make it > 10.

	product := 1.0;
	while (arg <= 10.0) do
	begin
		product := product * arg;
		arg := arg + 1.0;
	end;

// Use a polynomial approximation to Stirling's formula.
// N.B. The real Stirling's formula is used here, not the simpler, but less
//     accurate formula given by De Moivre in a letter to Stirling, which
//     is the one usually quoted.

	arg := arg - 0.5;
	temp := 1.0 / Sqr (arg);
	Result := lnrt2pi + arg * (Ln (arg) - 1.0 +
		(((a4 * temp + a3) * temp + a2) * temp + a1) * temp) - Ln (product);

	if reflect then
	begin
		temp := Sin (pi * x);
		Result := Ln (pi / temp) - Result;
	end;
end;

{     Calculate a binomial probability
}
function bin_prob (const n: Int64; const p: Extended; const r: Int64): Extended;
begin
	Result := Exp (lngamma (n +1.0) - lngamma (r + 1.0) - lngamma (n - r + 1.0)
		+ r * Ln (p) + (n - r) * Ln (1.0 - p));
end;


{  FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
  This algorithm is suitable when many random variates are required
  with the SAME parameter values for n & p.

	P = BERNOULLI SUCCESS PROBABILITY
		  (0 <= REAL <= 1)
	N = NUMBER OF BERNOULLI TRIALS
		  (1 <= INTEGER)

  Reference: Kemp, C.D. (1986). `A modal method for generating binomial
		   variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
}
function  Random_Binomial1 (const n: Integer; const p: Extended;
	RandomGenerator: TRandomGenFunction): Int64; 
var
	ru, rd: Int64;
	r0: Int64;
	u, pd, pu: Real;
	odds_ratio, p_r: Real;
begin
	r0 := Trunc ((n + 1) * p);
	p_r := bin_prob (n, p, r0);
	if p < 1 then
		odds_ratio := p / (1.0 - p)
	else
		odds_ratio := VLarge;

	u := RandomGenerator;
	u := u - p_r;
	if (u < 0.0) then
	begin
		Result := r0;
		Exit;
	end;

	pu := p_r;
	ru := r0;
	pd := p_r;
	rd := r0;

	repeat
		Dec (rd);
		if (rd >= 0) then
		begin
			pd := pd * (rd + 1.0) / (odds_ratio * (n - rd));
			u := u - pd;
			if (u < 0.0) then
			begin
				Result := rd;
				Exit;
			end;
		end;

		Inc (ru);
		if (ru <= n) then
		begin
			pu := pu * (n - ru+ 1.0) * odds_ratio / ru;
			u := u - pu;
			if (u < 0.0) then
			begin
				Result := ru;
				Exit;
			end;
		end;
	until False;
end;

function  Random_Binomial1 (const n: Integer; const p: Extended): Int64;
begin
	Result := Random_Binomial1 (n, p, DelphiRandom);
end;
{**********************************************************************
	 Translated to Fortran 90 by Alan Miller from:
						 RANLIB

	 Library of Fortran Routines for Random Number Generation

				   Compiled and Written by:

					   Barry W. Brown
					    James Lovato

			 Department of Biomathematics, Box 237
			 The University of Texas, M.D. Anderson Cancer Center
			 1515 Holcombe Boulevard
			 Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

				 GENerate BINomial random deviate

						 Function

	 Generates a single random deviate from a binomial
	 distribution whose number of trials is N and whose
	 probability of an event in each trial is P.

						 Arguments

	 N  --> The number of trials in the binomial distribution
		   from which a random deviate is to be generated.
						 INTEGER N

	 P  --> The probability of an event in each trial of the
		   binomial distribution from which a random deviate
		   is to be generated.
						 REAL P

!     FIRST --> Set FIRST = .TRUE. for the first call to perform initialization
!               the set FIRST = .FALSE. for further calls using the same pair
!               of parameter values (N, P).
!                              LOGICAL FIRST

!     random_binomial2 <-- A random deviate yielding the number of events
!                from N independent trials, each of which has
!                a probability of event P.
!                              INTEGER random_binomial

!                              Method

!     This is algorithm BTPE from:

!         Kachitvichyanukul, V. and Schmeiser, B. W.
!         Binomial Random Variate Generation.
!         Communications of the ACM, 31, 2 (February, 1988) 216.

!**********************************************************************
}
function Random_Binomial2 (const n: Int64; const pp: Extended;
	RandomGenerator: TRandomGenFunction): Int64; 
var
	alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2: Extended;
	i: Integer;
	ix, ix1, k, mp: Int64;
	m: Int64;
	p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,
	xlr, p2, p3, p4, qn, r, g: Extended;
	Done: Boolean;
begin
	// SETUP
	p := Min (pp, 1.0 - pp);
	q := 1.0 - p;
	xnp := n * p;

	if (xnp > 30.0) then
	begin
		ffm := xnp + p;
		m := Trunc (ffm);
		fm := m;
		xnpq := xnp * q;
		p1 := INT (2.195 * Sqrt (xnpq) - 4.6 * q) + 0.5;
		xm := fm + 0.5;
		xl := xm - p1;
		xr := xm + p1;
		c := 0.134 + 20.5 / (15.3 + fm);
		al := (ffm - xl) / (ffm - xl * p);
		xll := al * (1.0 + 0.5 * al);
		al := (xr - ffm) / (xr * q);
		xlr := al * (1.0 + 0.5 * al);
		p2 := p1 * (1.0 + c + c);
		p3 := p2 + c / xll;
		p4 := p3 + c / xlr;

		// GENERATE VARIATE, Binomial mean at least 30.

		ix := 0;
		repeat;
			u := RandomGenerator; //20
			u := u * p4;
			v := RandomGenerator;

			// TRIANGULAR REGION
			if (u <= p1) then
			begin
				ix := Trunc (xm - p1 * v + u);
				Break;
			end;

			// PARALLELOGRAM REGION
			if  (u <= p2) then
			begin
				x := xl + (u -p1) / c;
				v := v * c + 1.0 - Abs (xm - x) / p1;
				if (v > 1.0) or (v <= 0) then
					Continue;
				ix := Trunc (x);
			end
			else
			begin
				// LEFT TAIL
				if (u <= p3) then
				begin
					ix := Trunc (xl + Ln (v) / xll);
					if (ix < 0) then
						Continue;
					v := v * (u-p2) * xll
				end
				// RIGHT TAIL
				else
				begin
					ix := Trunc (xr - Ln (v) / xlr);
					if (ix > n) then
						Continue;
					v := v * (u-p3) * xlr;
				end;
			end;

			// DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST

			k := Abs (ix - m);
			if (k <= 20) or (k >= xnpq / 2 - 1) then
			begin
				// EXPLICIT EVALUATION
				f := 1.0;
				r := p / q;
				g := (n + 1) * r;
				if (m < ix) then
				begin
					mp := m + 1;
					for i := mp to ix do
						f := f * (g / i - r);
				end
				else if (m > ix) then
				begin
					ix1 := ix + 1;
					for i := ix1 to m do
						f := f / (g / i - r);
				end;

				if (v > f) then
					Continue
				else
					Break
			end;

			// SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))

			amaxp := (k / xnpq) * ((k * (k / 3.0 + 0.625)
				+ 0.1666666666666) / xnpq + 0.5);
			ynorm := - Sqr (k) / (2.0 * xnpq);
			alv := Ln (v);
			if (alv < ynorm - amaxp) then
				Break;
			if (alv > ynorm + amaxp) then
				Continue;
			// STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
			//  THE FINAL ACCEPTANCE/REJECTION TEST

			x1 := ix + 1;
			f1 := fm + 1.0;
			z := n + 1 - fm;
			w := n - ix + 1.0;
			z2 := Sqr (z);
			x2 := Sqr (x1);
			f2 := Sqr (f1);
			w2 := Sqr (w);
			if (alv - (xm * Ln (f1 / x1) + (n - m + 0.5) * Ln (z / w)
				+ (ix - m) * Ln (w * p / (x1 * q))
				+ (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0
				+ (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0
				+ (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0
				+ (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.0)
				> 0.0) then
			begin
				Continue;
			end
			else
				Break;
		until False;
	end
	else
	begin
		// INVERSE CDF LOGIC FOR MEAN LESS THAN 30
		qn := Power (q, n);
		r := p / q;
		g := r * (n + 1);

		Done := False;
		repeat
			ix := 0;
			f := qn;
			u := RandomGenerator;
			while (u >= f) do
			begin
				if (ix > 110) then
					Done := False
				else
				begin
					Done := True;
					u := u - f;
					ix := ix + 1;
					f := f * (g / ix - r);
				end;
			end;
		until Done;
	end;

	if (pp > 0.5) then
		ix := n - ix;

	Result := ix;
end;

function Random_Binomial2 (const n: Int64; const pp: Extended): Int64;
begin
	Result := Random_Binomial2 (n, pp);
end;

{ Adapted from Fortran 77 code from the book:
	 Dagpunar, J. 'Principles of random variate generation'
	 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
  INVERSION AND/OR THE REPRODUCTIVE PROPERTY.

	SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
	   = the `power' parameter of the negative binomial
		  (0 < REAL)
	P = BERNOULLI SUCCESS PROBABILITY
		  (0 < REAL < 1)

  THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
  OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
  THE REPRODUCTIVE PROPERTY IS USED.
}
function Random_Neg_Binomial (const sk, p: Extended;
	RandomGenerator: TRandomGenFunction): Int64;
const
	h = 0.7;
var
	q, x, st, uln, v, r, s, y, g: Extended;
	k, n: Int64;
	i: Integer;
begin
	if  (sk <= 0.0) or (p <= 0.0) or (p >= 1.0) then
		raise EMathError.Create (rsInvalidValue);

	q := 1.0 - p;
	x := 0.0;
	st := sk;
	if  (p > h) then
	begin
		v := 1.0 / Ln (p);
		k := Trunc (st);
		for i := 1 to k do
		begin
			repeat
				r := RandomGenerator;
			until r > 0.0;
			n := Trunc (v * Ln (r));
			x := x + n;
		end;
		st := st - k;
	end;

	s := 0.0;
	uln := -Ln (VSmall);
	if (st > -uln / Ln (q)) then
		raise EMathError.Create (rsInvalidValue);

	y := Power (q, st);
	g := st;
	r := RandomGenerator;
	while (y <= r) do
	begin
		r := r - y;
		s := s + 1.0;
		y := y * p * g / s;
		g := g + 1.0;
	end;

	Result := Trunc (x + s + 0.5)
end;

function Random_Neg_Binomial (const sk, p: Extended): Int64;
begin
	Result := Random_Neg_Binomial (sk, p, DelphiRandom);
end;

//     Gaussian integration of exp(k.cosx) from a to b.
	const dk: Extended);
var
	xmid, range, x1, x2: Extended;
	X, W: array [1..3] of Extended;
	I: Integer;
begin
	X [1] := 0.238619186083197;
	X [2] := 0.661209386466265;
	X [3] := 0.932469514203152;
	W [1] := 0.467913934572691;
	W [2] := 0.360761573048139;
	W [3] := 0.171324492379170;

	xmid := (a + b) / 2.0;
	range := (b - a) / 2.0;

	res := 0.0;
	for I := 1 to 3 do
	begin
		x1 := xmid + X  [I] * range;
		x2 := xmid - X [I] * range;
		res := res + W [I] * (Exp (dk * Cos (x1)) + Exp (dk * Cos (x2)))
	end;

	res := res * range;
end;

{     Algorithm VMD from:
	 Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
	 comparison of random numbers', J. of Appl. Statist., 17, 165-168.

	 Fortran 90 code by Alan Miller
	 CSIRO Division of Mathematical & Information Sciences

	 Arguments:
	 k (real)        parameter of the von Mises distribution.
}

function Random_von_Mises (const k: Extended;
	RandomGenerator: TRandomGenFunction): Extended;
var
	j, n, jj: Integer;
	nk: Integer;
	p: array [1..20] of Extended;
	theta: array [0..20] of Extended;
	xx, sump, r, th, lambda, rlast: Extended;
	dk: Extended;
begin
	if (k < 0.0) then
		raise EMathError.Create (rsInvalidValue);

	nk := Trunc (k + k + 1.0);
	if (nk > 20) then
		raise EMathError.Create (rsInvalidValue);

	dk := k;
	theta [0] := 0.0;
	if (k > 0.5) then
	begin
		// Set up array p of probabilities.
		sump := 0.0;
		for j := 1 to nk do
		begin
			if (j < nk) then
				theta [j] := ArcCos (1.0 - j / k)
			else
				theta [nk] := Pi;

			// Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
			Integral (theta [j - 1], theta [j], p[j], dk);
			sump := sump + p [j];
		end;
		for j := 1 to nk do
		    p [j] := p [j] / sump;
	end
	else
	begin
		p [1] := 1.0;
		theta [1] := Pi;
	end;

	r := RandomGenerator;
	jj := nk;
	for j := 1 to nk do
	begin
		r := r - p [j];
		if (r < 0.0) then
		begin
			jj := j;
			Break;
		end;
	end;
	r := -r / p[jj];

	repeat
		th := theta [jj - 1] + r * (theta [jj] - theta [j - 1]);
		lambda := k - jj + 1.0 - k * Cos (th);
		n := 1;
		rlast := lambda;

		repeat
			r := RandomGenerator;
			if r <= rlast then
			begin
				Inc (n);
				rlast := r;
			end;
		until (r > rlast);

		if not Odd (n) then
			r := RandomGenerator
	until Odd (n);

	th := Abs (th);
	xx := (r - rlast) / (1.0 - rlast) - 0.5;
	if xx < 0 then
		Result := -1 * th
	else
		Result := th;
end;

function Random_von_Mises (const k: Extended): Extended;
begin
	Result := Random_von_Mises (k, DelphiRandom);
end;

{      Generate a random deviate from the standard Cauchy distribution
}
function Random_Cauchy (RandomGenerator: TRandomGenFunction): Extended;
var
	V1, V2: Extended;
begin
	repeat
		V1 := 2.0 * (RandomGenerator - 0.5);
		V2 := 2.0 * (RandomGenerator - 0.5);
	until (Abs (V2) > VSmall) and (Sqr (V1) + Sqr (V2) < 1.0);
	Result := V1 / V2;
end;

function Random_Cauchy: Extended;
begin
	Result := Random_Cauchy (DelphiRandom);
end;

{**********************************************************************
	 Translated to Fortran 90 by Alan Miller from:
					   RANLIB

	 Library of Fortran Routines for Random Number Generation

				 Compiled and Written by:

					 Barry W. Brown
					  James Lovato

		    Department of Biomathematics, Box 237
		    The University of Texas, M.D. Anderson Cancer Center
		    1515 Holcombe Boulevard
		    Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

				 GENerate POIsson random deviate

					    Function

  Generates a single random deviate from a Poisson distribution with mean mu.

					    Arguments

	 mu --> The mean of the Poisson distribution from which
		   a random deviate is to be generated.
						 REAL mu

						 Method

	 For details see:

			 Ahrens, J.H. and Dieter, U.
			 Computer Generation of Poisson Deviates
			 From Modified Normal Distributions.
			 ACM Trans. Math. Software, 8, 2
			 (June 1982),163-179

	 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
	 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
}
function Random_Poisson(var mu: Extended;
	RandomGenerator: TRandomGenFunction): Int64;
const
	a0 = -0.5;
	a1 = 0.3333333;
	a2 = -0.2500068;
	a3 = 0.2000118;
	a4 = -0.1661269;
	a5 = 0.1421878;
	a6 = -0.1384794;
	a7 = 0.1250060;

var
	b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,
	omega, px, py, t, u, v, x, xx: Extended;
	s, d, p, q, p0: Extended;
	j, k, kflag: Integer;
	l, m: Integer;
	pp: array [1..35] of Extended;
	fact: array [1..10] of Extended;
	FirstK0: Boolean;
begin
	u := 0;
	e := 0;
	fk := 0;
	difmuk := 0;
	
	fact [1] := 1; // Factorial 0
	fact [2] := 1;
	fact [3] := 2;
	fact [4] := 6;
	fact [5] := 24;
	fact [6] := 120;
	fact [7] := 720;
	fact [8] := 5040;
	fact [9] := 40320;
	fact [10] := 362880; // Factorial 9

	Result := 0;

	if (mu > 10.0) then
	begin
		// C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)

		FirstK0 := False;

		s := Sqrt (mu);
		d := 6.0 * Sqr (mu);

		// THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
		// PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
		//   IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

		l := Trunc (mu - 1.1484);

		// STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

		g := mu + s * Random_Normal (RandomGenerator);
		if (g > 0.0) then
		begin
			Result := Trunc (g);

			// STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH

			if (Result >= l) then
				Exit;

			// STEP S. SQUEEZE ACCEPTANCE - SAMPLE U

			fk := Result;
			difmuk := mu - fk;
			u := RandomGenerator;
			if (d * u >= Sqr (difmuk) * difmuk) then
				Exit;
		end;

		// STEP P. PREPARATIONS FOR STEPS Q AND H.
		// (RECALCULATIONS OF PARAMETERS IF NECESSARY)
		// .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
		// THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
		// APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
		// C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

		omega := 0.3989423 / s;
		b1 := 0.4166667E-1 / mu;
		b2 := 0.3 * Sqr (b1);
		c3 := 0.1428571 * b1 * b2;
		c2 := b2 - 15.0 * c3;
		c1 := b1 - 6.0 * b2 + 45.0 * c3;
		c0 := 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
		c := 0.1069 / mu;

		kflag := -1;

		if (g >= 0.0) then
		begin

			// 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
			kflag := 0;
			FirstK0 := True;
		end;

		// STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
		// DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
		// (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

		repeat // 50
			if not FirstK0 then
			begin
				repeat // 50
					e := Random_Exponential (RandomGenerator);
					u := RandomGenerator;
					u := u + u - 1.0;
					if u < 0 then
						t := 1.8 - abs (e)
					else
						t := 1.8 + abs (e);
				until t > -0.6744;

				Result := Trunc (mu + s * t);
				fk := Result;
				difmuk := mu - fk;

				// 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

				kflag := 1
			end
			else
				FirstK0 := False;

			// STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
			// CASE ival < 10 USES FACTORIALS FROM TABLE FACT

			if (Result < 10) then // 70
			begin
				px := -mu;
				py := Power (mu, Result) / fact [Result + 1];
			end
			else
			begin
				// CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
				//  A0-A7 FOR ACCURACY WHEN ADVISABLE
				//  .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

				del := 0.8333333E-1 / fk; // 80
				del := del - 4.8 * Sqr (del) * del;
				v := difmuk / fk;
				if Abs (v) > 0.25 then
					px := fk * Ln (1.0 + v) - difmuk - del
				else
					px := fk * Sqr (v) * (((((((a7 * v + a6) * v + a5) * v + a4)
						* v + a3) * v + a2) * v + a1) * v + a0) - del;
				  py := 0.3989423 / Sqrt (fk);
			end;
			x := (0.5 - difmuk) / s; // 110
			xx := Sqr (x);
			fx := -0.5 * xx;
			fy := omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
			if (kflag <= 0) then
			begin
				// STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

				if (fy - u * fy <= py *Exp (px - fx)) then // 40
					Exit;
			end
			else
			begin
			// STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

				if (c * Abs (u) <= py * Exp (px + e) - fy * Exp (fx + e)) then // 60
					Exit;
			end;
		until False;
	end
	else
	begin

		// C A S E  B.    mu < 10
		// START NEW TABLE AND CALCULATE P0 IF NECESSARY

		m := MAX (1, Trunc (mu));
		l := 0;
		p := Exp (-mu);
		q := p;
		p0 := p;

		//  STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

		repeat
			u := RandomGenerator;
			Result := 0;
			if (u <= p0) then
				Exit;

			// STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
			// PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
			// (0.458=PP(9) FOR MU=10)

			if  l <> 0 then
			begin
				j := 1;
				if (u > 0.458) then
					 j := MIN (l, m);
				for k := j to l do
					if (u <= pp [k]) then
					begin
						Result := K;
						Exit;
					end;

				if l = 35 then
					Continue;
			end;

			// STEP C. CREATION OF NEW POISSON PROBABILITIES P
			// AND THEIR CUMULATIVES Q=PP(K)

			l := l + 1; // 150
			for k := l to 35 do
			begin
				p := p * mu / k;
				q := q + p;
				pp [k] := q;
				if (u <= q) then
				begin
					Result := K;
					Exit;
				end;
			end;
			l := 35;
		until False
	end;
end;

function Random_Poisson(var mu: Extended): Int64;
begin
	Result := Random_Poisson (mu, DelphiRandom);
end;


end.


.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]