| Random_Binomial1 Routines |
Unit
QESBPCSRandom
| Overloaded Variants |
| Function Random_Binomial1(const n: Integer; const p: Extended): Int64; |
| Function Random_Binomial1(const n: Integer; const p: Extended; RandomGenerator: TRandomGenFunction): Int64; |
Declaration
Function Random_Binomial1(const n: Integer; const p: Extended): Int64;
Description
This algorithm is suitable when many random variates are required with the SAME parameter values for n & p. Reference: Kemp, C.D. (1986). `A modal method for generating binomial variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
| Parameters |
| n | Number of Bernoulli Trials, must be positive. |
| p | Bernoulli Probability of Success must be in [0, 1]. |
| 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_Binomial1 (const n: Integer; const p: Extended): Int64;
begin
Result := Random_Binomial1 (n, p, DelphiRandom);
End; |
Declaration
Function Random_Binomial1(const n: Integer; const p: Extended; RandomGenerator: TRandomGenFunction): Int64;Implementation
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; |
|
|