stblinv.m

matlab alpha稳定分布

function x = stblinv(u,alpha,beta,gam,delta,varargin)
%X = STBLINV(U,ALPHA,BETA,GAM,DELTA) returns values of the inverse CDF of
% the alpha-stable distribution with characteristic exponent ALPHA, skewness
% BETA, scale GAM, and location DELTA at the values in the array U.   
%
% This alogrithm uses a combination of Newton's method and the bisection 
% method to compute the inverse cdf to a tolerance of 1e-6;
%
% X = STBLINV(U,ALPHA,BETA,GAM,DELTA,'quick') returns a linear interpolated
% approximation of the inverse CDF based on a table of pre-calculated
% values.  The table contains exact values at
%   ALPHA = [.1 : .1: 2]
%   BETA  = [0: .2 :  1]
%   U     = [ .1 : .1 : .9 ]
% If U < .1 or U > .9, the 'quick' option approximates the CDF with its 
% asymptotic form which is given in [1], page 16, Property 1.2.15.  Results
% for U outside the interval [.1:.9] may vary.
% 
% If abs(ALPHA - 1) < 1e-5,  ALPHA is rounded to 1.
%
% See also:  STBLRND, STBLPDF, STBLCDF, STBLFIT
%
% Reference:
%  [1] G. Samorodnitsky & M. S. Taqqu (1994)
%      "Stable Non-Gaussian Random Processes, Stochastic Models with 
%       Infinite Variance"  Chapman & Hall
%
if nargin < 5
    error('stblcdf:TooFewInputs','Requires at least five input arguments.'); 
end

% Check parameters
if alpha <= 0 || alpha > 2 || ~isscalar(alpha)
    error('stblcdf:BadInputs',' "alpha" must be a scalar which lies in the interval (0,2]');
end
if abs(beta) > 1 || ~isscalar(beta)
    error('stblcdf:BadInputs',' "beta" must be a scalar which lies in the interval [-1,1]');
end
if gam < 0 || ~isscalar(gam)
    error('stblcdf:BadInputs',' "gam" must be a non-negative scalar');
end
if ~isscalar(delta)
    error('stblcdf:BadInputs',' "delta" must be a scalar');
end

if (1e-5 < abs(alpha - 1)  && abs(alpha - 1) < .02) || alpha < .02
    warning('stblcdf:ScaryAlpha',...
        'Difficult to approximate cdf for alpha close to 0 or 1')
end

quick = false;
if ~isempty(varargin)
    if strcmp(varargin{1},'quick')
        quick = true;
    end
end

% For Newton's Method
itermax = 30;
maxbisecs = 30;
tol = 1e-8;

% Return NaN for out of range parameters or probabilities.
u(u < 0 | 1 < u) = NaN;

% Check to see if you are in a simple case, if so be quick.
if alpha == 2                  % Gaussian distribution 
    x = - 2 .* erfcinv(2*u);
    x = x*gam + delta;    %  
elseif alpha==1 && beta == 0   % Cauchy distribution
    x =  tan(pi*(u - .5) );   
    x = x*gam + delta;
elseif alpha == .5 && abs(beta) == 1 % Levy distribution 
    x = .5 * beta ./(erfcinv(u)).^2;
    x = x*gam + delta;
else                % Gen. Case
    
    % Flip sign of beta if necessary
    if beta < 0
        signBeta = -1;
        u = 1-u;
        beta = -beta;
    else
        signBeta = 1;
    end
    
    % Calculate additional shift for (M) parameterization 
    if abs(alpha - 1) > 1e-5
        deltaM = -beta * tan(alpha*pi/2); % 
    else
        deltaM = 0;
    end
     
    x = intGuess(alpha,beta,u);
    if ~quick
        % Newton's Method
        F = stblcdf(x,alpha,beta,1,deltaM) - u;
        diff = max(abs(F),0); % max incase of NaNs
        bad = diff > tol;
        iter = 1;
        Fold = F;
        xold = x; 
        while any(bad) && iter < itermax 
            
            % Perform Newton step 
            % If Fprime = 0, step closer to origin instead
            Fprime = stblpdf(x(bad),alpha,beta,1,deltaM,1e-8);
            x(bad) = x(bad) - F(bad) ./ Fprime;
            blowup = isinf(x) | isnan(x);
            if ~isempty(blowup) 
                x(blowup) = xold(blowup) / 2;
            end

            F(bad) = stblcdf(x(bad),alpha,beta,1,deltaM) - u(bad);

            % Make sure we are getting closer, if not, do bisections until
            % we do.
            nocvg = abs(F) > 1.1*abs(Fold);
            bisecs = 0;
            while any(nocvg) && (bisecs < maxbisecs)
                x(nocvg) = .5*(x(nocvg) + xold(nocvg));
                F(nocvg) = stblcdf(x(nocvg),alpha,beta,1,deltaM) - u(nocvg);
                nocvg = abs(F) > 1.1*abs(Fold);
                bisecs = bisecs + 1;
            end

            % Test for convergence
            diff = max(abs(F),0); % max incase of NaNs
            bad = diff > tol;

            % Save for next iteration
            xold = x;
            Fold = F;
            iter = iter + 1;
        end
    end
      
    % Un-standardize
    if abs(1 - alpha) > 1e-5
        x = signBeta*(x - deltaM)*gam + delta;
    else
        x = signBeta*(x*gam + (2/pi) * beta * gam * log(gam)) + delta;
    end
       
end


end



%===================================================================
%===================================================================
%===================================================================

function X0 = intGuess(alpha,beta,u)
% Look-up table of percentiles of standard stable distributions
% If .1 < u < .9, Interpolatates tabulated values to obtain initial guess
% If u < .1 or u > .9 uses asymptotic formulas to make a starting guess