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