Random_Inv_Gauss Routines |
Unit
QESBPCSRandom
Overloaded Variants |
Function Random_Inv_Gauss(const h, b: Extended): Extended; |
Function Random_Inv_Gauss(const h, b: Extended; RandomGenerator: TRandomGenFunction): Extended; |
Declaration
Function Random_Inv_Gauss(const h, b: Extended): Extended;
Description
Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9
Parameters |
h | Parameter of Distribution, must be positive. |
b | Parameter of Distribution, must be positive. |
RandomGenerator | Optional Function to use for Uniform Random Number Generator. If omitted, Delphi's Random function is used, and if this is done remember to call Randomize if you don't want repeated values. |
Category
Arithmetic Routines for FloatsImplementation
function Random_Inv_Gauss (const h, b: Extended): Extended; begin Result := Random_Inv_Gauss (h, b, DelphiRandom); End; |
Declaration
Function Random_Inv_Gauss(const h, b: Extended; RandomGenerator: TRandomGenFunction): Extended;Implementation
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 := XtoY (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; |
|