Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

vonmises_cdf(x, kappa, mu)
% vonmises_cdf.m - evaluates a Cumulative Von MisesDistribution.
%   See "ACM TOMS Alg 518", ACM Trans Math SW, v3, n3, Sept 1977, p. 279.
%
%  Created by Jim Huntley,  11/06/03
%

function [cdf] = vonmises_cdf(x, kappa, mu)

% Initializations.
a1 = 12;
a2 = 0.8;
a3 = 8;
a4 = 1;
ck = 10.5;
c1 = 56;


if(x-mu <= -pi)
    cdf = 0;
    return
end
if(x(jx)-mu >= pi)
    cdf(jx) = 1;
    return
end

% Convert x-mu to the range 0 -> 2*pi.
z = kappa;
u = mod(x-mu+pi,2*pi);
if(u < 0)
    u = u + 2*pi;
end
y = u - pi;

% For small values of kappa,
if(z <= ck)
    v = 0;
    if(z > 0)
        ip = fix(z*a2 - a3/(z+a4) + a1);
        p = ip;
        %ip = fix(p);
        sy = sin(y);
        cy = cos(y);
        y = p * y;
        sny = sin(y);
        cny = cos(y);
        r = 0;
        z = 2 / z;
        for jn = 2:ip
            p = p - 1;
            y = sny;
            sny = sny*cy - cny*sy;
            cny = cny*cy + y*sy;
            r = 1 / (p*z + r);
            v = (sny/p + v) * r;
        end
    end

    cdf = (0.5*u + v) / pi;
% For large values of kappa,   
else 

    c = 24 * z;
    v = c - c1;
    r = sqrt((54/(347/v + 26 - c) -6 + c)/12);
    z = sin(0.5*y) * r;
    s = 2*z^2;
    v = v - s + 3;
    y = (c - s - s -16)/3;
    y = ((s + 1.75)*s + 83.5)/v - y;
    arg = z*(1 - s/y^2);
    erfx = erf(arg);
    cdf = 0.5*erfx + 0.5;
end

cdf = max(cdf,0);
cdf = min(cdf,1);


return

Contact us