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