Random_Beta Routines |
Unit
QESBPCSRandom
Overloaded Variants |
Function Random_Beta(const aa, bb: Extended): Extended; |
Function Random_Beta(const aa, bb: Extended; RandomGenerator: TRandomGenFunction): Extended; |
Declaration
Function Random_Beta(const aa, bb: Extended): Extended;
Description
Uses Cheng'S log logistic method. 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). Uses Cheng'S log logistic method.
Parameters |
aa | Shape Parameter for Beta Distribution - must be positive. |
bb | Shape Parameter for Beta 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_Beta (const aa, bb: Extended): Extended; begin Result := Random_Beta (aa, bb, DelphiRandom); End; |
Declaration
Function Random_Beta(const aa, bb: Extended; RandomGenerator: TRandomGenFunction): Extended;Implementation
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 + XtoY (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 > XtoY (1.0 + 1.0 / d, f) then Continue; Result := 1.0 end; Done := True until Done; if Swap then Result := 1 - Result; End; |
|