Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

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