| Random_von_Mises Routines |
Unit
QESBPCSRandom
| Overloaded Variants |
| Function Random_von_Mises(const k: Extended): Extended; |
| Function Random_von_Mises(const k: Extended; RandomGenerator: TRandomGenFunction): Extended; |
Declaration
Function Random_von_Mises(const k: Extended): Extended;
Description
Algorithm VMD from: Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a comparison of random numbers', J. of Appl. Statist., 17, 165-168.
Fortran 90 code by Alan Miller
CSIRO Division of Mathematical & Information Sciences
| Parameters |
| k | Parameter of the von Mises distribution. |
| 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_von_Mises (const k: Extended): Extended;
begin
Result := Random_von_Mises (k, DelphiRandom);
End; |
Declaration
Function Random_von_Mises(const k: Extended; RandomGenerator: TRandomGenFunction): Extended;Implementation
function Random_von_Mises (const k: Extended;
RandomGenerator: TRandomGenFunction): Extended;
var
j, n, jj: Integer;
nk: Integer;
p: array [1..20] of Extended;
theta: array [0..20] of Extended;
xx, sump, r, th, lambda, rlast: Extended;
dk: Extended;
begin
if (k < 0.0) then
raise EMathError.Create (rsInvalidValue);
nk := Trunc (k + k + 1.0);
if (nk > 20) then
raise EMathError.Create (rsInvalidValue);
dk := k;
theta [0] := 0.0;
if (k > 0.5) then
begin
// Set up array p of probabilities.
sump := 0.0;
for j := 1 to nk do
begin
if (j < nk) then
theta [j] := ESBArcCos (1.0 - j / k)
else
theta [nk] := ESBPi;
// Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
Integral (theta [j - 1], theta [j], p [j], dk);
sump := sump + p [j];
end;
for j := 1 to nk do
p [j] := p [j] / sump;
end
else
begin
p [1] := 1.0;
theta [1] := ESBPi;
end;
r := RandomGenerator;
jj := nk;
for j := 1 to nk do
begin
r := r - p [j];
if (r < 0.0) then
begin
jj := j;
Break;
end;
end;
r := -r / p [jj];
repeat
th := theta [jj - 1] + r * (theta [jj] - theta [j - 1]);
lambda := k - jj + 1.0 - k * Cos (th);
n := 1;
rlast := lambda;
repeat
r := RandomGenerator;
if r <= rlast then
begin
Inc (n);
rlast := r;
end;
until (r > rlast);
if not Odd (n) then
r := RandomGenerator
until Odd (n);
th := Abs (th);
xx := (r - rlast) / (1.0 - rlast) - 0.5;
if xx < 0 then
Result := -1 * th
else
Result := th;
End; |
|
|