| 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; |
|
|