Random_von_Mises Routines
Function generates a Random von Mises Variate.

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

Implementation

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;


HTML generated by Time2HELP
http://www.time2help.com