stabrnd.m

A toolbox for the non-local means algorithm

function [x] = stabrnd(alpha, beta, c, delta, m, n)

% STABRND.M
% Stable Random Number Generator (McCulloch 12/18/96)
%
%   x = stabrnd(alpha, beta, c, delta, m, n);
%
% Returns m x n matrix of iid stable random numbers with 
%   characteristic exponent alpha in [.1,2], skewness parameter
%   beta in [-1,1], scale c > 0, and location parameter delta.
% Based on the method of J.M. Chambers, C.L. Mallows and B.W.
%   Stuck, "A Method for Simulating Stable Random Variables," 
%   JASA 71 (1976): 340-4.  
% Encoded in MATLAB by J. Huston McCulloch, Ohio State
%   University Econ. Dept. ([e-mail]).  This 12/18/96
%   version uses 2*m*n calls to RAND, and does not rely on 
%   the STATISTICS toolbox.
% The CMS method is applied in such a way that x will have the 
%   log characteristic function 
%        log E exp(ixt) = i*delta*t + psi(c*t), 
%   where
%     psi(t) = -abs(t)^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))
%                              for alpha ~= 1,
%            = -abs(t)*(1+i*beta*(2/pi)*sign(t)*log(abs(t))),
%                              for alpha = 1.
% With this parameterization, the stable cdf S(x; alpha, beta, 
%   c, delta) equals S((x-delta)/c; alpha, beta, 1, 0).  See my 
%   "On the parametrization of the afocal stable distributions,"
%   _Bull. London Math. Soc._ 28 (1996): 651-55, for details.
% When alpha = 2, the distribution is Gaussian with mean delta 
%   and variance 2*c^2, and beta has no effect.
% When alpha > 1, the mean is delta for all beta.  When alpha 
%   <= 1, the mean is undefined.
% When beta = 0, the distribution is symmetrical and delta is 
%   the median for all alpha.  When alpha = 1 and beta = 0, the 
%   distribution is Cauchy (arctangent) with median delta.
% When the submitted alpha is > 2 or < .1, or beta is outside 
%   [-1,1], an error message is generated and x is returned as a 
%   matrix of NaNs.
% Alpha < .1 is not allowed here because of the non-negligible 
%   probability of overflows.  
%
% If you're only interested in the symmetric cases, you may just 
%   set beta = 0 and skip the following considerations:
% When beta > 0 (< 0), the distribution is skewed to the right 
%   (left).
% When alpha < 1, delta, as defined above, is the unique fractile
%   that is invariant under averaging of iid contributions.  I 
%   call such a fractile a "focus of stability."  This, like the 
%   mean, is a natural location parameter.
% When alpha = 1, either every fractile is a focus of stability, 
%   as in the beta = 0 Cauchy case, or else there is no focus of 
%   stability at all, as is the case for beta ~=0.  In the latter
%   cases, which I call "afocal," delta is just an arbitrary 
%   fractile that has a simple relation to the c.f.
% When alpha > 1 and beta > 0, med(x) must lie very far below 
%   the mean as alpha approaches 1 from above.  Furthermore, as 
%   alpha approaches 1 from below, med(x) must lie very far above
%   the focus of stability when beta > 0.  If beta ~= 0, there  
%   is therefore a discontinuity in the distribution as a function
%   of alpha as alpha passes 1, when delta is held constant.
% CMS, following an insight of Vladimir Zolotarev, remove this
%   discontinuity by subtracting 
%          beta*c*tan(pi*alpha/2)
%   (equivalent to their -tan(alpha*phi0)) from x for alpha ~=1
%   in their program RSTAB, a.k.a. RNSTA in IMSL (formerly GGSTA).
%   The result is a random number whose distribution is a contin-
%   uous function of alpha, but whose location parameter (which I 
%   call zeta) is a shifted version of delta that has no known 
%   interpretation other than computational convenience.  
%   The present program restores the more meaningful "delta"  
%   parameterization by using the CMS (4.1), but with 
%   beta*c*tan(pi*alpha/2) added back in (ie with their initial
%   tan(alpha*phi0) deleted).  RNSTA therefore gives different 
%   results than the present program when beta ~= 0.  However,
%   the present beta is equivalent to the CMS beta' (BPRIME).
% Rather than using the CMS D2 and exp2 functions to compensate
%   for the ill-condition of the CMS (4.1) when alpha is very 
%   near 1, the present program merely fudges these cases by 
%   computing x from their (2.4) and adjusting for 
%   beta*c*tan(pi*alpha/2) when alpha is within 1.e-8 of 1. 
%   This should make no difference for simulation results with 
%   samples of size less than approximately 10^8, and then 
%   only when the desired alpha is within 1.e-8 of 1, but not 
%   equal to 1.
% The frequently used Gaussian and symmetric cases are coded 
%   separately so as to speed up execution.
%
% Additional references:
% V.M. Zolotarev, _One Dimensional Stable Laws_, Amer. Math. 
%   Soc., 1986.
% G. Samorodnitsky and M.S. Taqqu, _Stable Non-Gaussian Random
%   Processes_, Chapman & Hill, 1994.
% A. Janicki and A. Weron, _Simulaton and Chaotic Behavior of 
%   Alpha-Stable Stochastic Processes_, Dekker, 1994.
% J.H. McCulloch, "Financial Applications of Stable Distributons,"
%   _Handbook of Statistics_ Vol. 14, forthcoming early 1997.

% Errortraps:
if alpha < .1 | alpha > 2
    disp('Alpha must be in [.1,2] for function STABRND.')
    alpha
    x = NaN * zeros(m,n);
    return
end
if abs(beta) > 1
    disp('Beta must be in [-1,1] for function STABRND.')
    beta
    x = NaN * zeros(m,n);
    return
end

% Generate exponential w and uniform phi:
w = -log(rand(m,n));
phi = (rand(m,n)-.5)*pi;

% Gaussian case (Box-Muller):
if alpha == 2
    x = (2*sqrt(w) .* sin(phi));
    x = delta + c*x;
    return
end

% Symmetrical cases:
if beta == 0
    if alpha == 1   % Cauchy case
        x = tan(phi);
    else
        x = ((cos((1-alpha)*phi) ./ w) .^ (1/alpha - 1)    ...
            .* sin(alpha * phi) ./ cos(phi) .^ (1/alpha));
    end

    % General cases:
else
    cosphi = cos(phi);
    if abs(alpha-1) > 1.e-8
        zeta = beta * tan(pi*alpha/2);
        aphi = alpha * phi;
        a1phi = (1 - alpha) * phi;
        x = ((sin(aphi) + zeta * cos(aphi)) ./ cosphi)  ...
            .* ((cos(a1phi) + zeta * sin(a1phi))        ...
            ./ (w .* cosphi)) .^ ((1-alpha)/alpha);
    else
        bphi = (pi/2) + beta * phi;
        x = (2/pi) * (bphi .* tan(phi) - beta * log((pi/2) * w ...
            .* cosphi ./ bphi));
        if alpha ~= 1
            x = x + beta * tan(pi * alpha/2);
        end
    end
end

% Finale:
x = delta + c * x;
return
% End of STABRND.M