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