Code covered by the BSD License  

Highlights from
STBL: Alpha stable distributions for MATLAB

image thumbnail

STBL: Alpha stable distributions for MATLAB

by

 

Functions to handle one-dimensional alpha stable distributions

stblinv(u,alpha,beta,gam,delta,varargin)
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 

utemp = u(:);
X0 = zeros(size(utemp));
alpha = max(alpha,.1);
if beta == 1
    utemp(utemp < .1) = .1;  % bring these into table range 
end                          % since asyp. formulas don't apply if beta=1.
 
high = (utemp > .9);  
low = (utemp < .1);
middle = ~high & ~low;

% Use asymptotic formulas to guess high and low
if any(high | low)
    if alpha ~= 1 
        C = (1-alpha) / ( gamma(2 - alpha) * cos(pi*alpha/2) );
    else
        C = 2/pi;
    end
    X0(high) = ( (1-u(high))/(C * .5 * (1 + beta)) ).^(-1/alpha); 
    X0(low)  = -(u(low)/(C * .5 * (1 - beta))).^(-1/alpha);
end

% Use pre-calculated lookup table
if any(middle)
    [Alp Bet P] = meshgrid(.1:.1:2 , 0:.2:1 , .1:.1:.9 ); 
    stblfrac = zeros(6,20,9);
    stblfrac(:,1:5,1) = ...  % 
      [-1.890857122067030e+006   -1.074884919696010e+003   -9.039223076384694e+001   -2.645987890965098e+001   -1.274134564492298e+001;...
       -1.476366405440763e+005   -2.961237538429159e+002   -3.771873580263473e+001   -1.357404219788403e+001   -7.411052003232824e+000;...
       -4.686998894118387e+003   -5.145071882481552e+001   -1.151718246460839e+001   -5.524535336243413e+000   -3.611648531595958e+000;...
       -2.104710824345458e+001   -3.379418096823576e+000   -1.919928049616870e+000   -1.508399002681057e+000   -1.348510542803496e+000;...
       -1.267075422596289e-001   -2.597188113311268e-001   -4.004811495862077e-001   -5.385024279816432e-001   -6.642916520777534e-001;...
       -1.582153175255304e-001   -3.110425775503970e-001   -4.383733961816599e-001   -5.421475800719634e-001   -6.303884905318050e-001];
    stblfrac(:,6:10,1) = ...
      [-7.864009406553024e+000   -5.591791397752695e+000   -4.343949435866958e+000   -3.580521076832391e+000   -3.077683537175253e+000;...
       -4.988799898398770e+000   -3.787942909197120e+000   -3.103035515608863e+000   -2.675942594722292e+000   -2.394177022026705e+000;...
       -2.762379160216148e+000   -2.313577186902494e+000   -2.052416861482463e+000   -1.893403771865641e+000   -1.796585983161395e+000;...
       -1.284465355994317e+000   -1.267907903071982e+000   -1.279742001004255e+000   -1.309886183701422e+000   -1.349392554642457e+000;...
       -7.754208907962602e-001   -8.732998811318613e-001   -9.604322013853581e-001   -1.039287445657237e+000   -1.111986321525904e+000;...
       -7.089178961038225e-001   -7.814055112235459e-001   -8.502117698317242e-001   -9.169548634355569e-001   -9.828374636178471e-001];
    stblfrac(:,11:15,1) = ...
      [-2.729262880847457e+000   -2.479627528870857e+000   -2.297138304998905e+000   -2.162196365947914e+000   -2.061462692277420e+000;...
       -2.202290611202918e+000   -2.070075681428623e+000   -1.979193969170630e+000   -1.917168989568703e+000   -1.875099179801364e+000;...
       -1.740583121589162e+000   -1.711775396141753e+000   -1.700465158047576e+000   -1.700212465596452e+000   -1.707238269631509e+000;...
       -1.391753942957071e+000   -1.434304119387730e+000   -1.476453646904256e+000   -1.518446568503842e+000   -1.560864595722380e+000;...
       -1.180285915835185e+000   -1.245653509438976e+000   -1.309356535558631e+000   -1.372547245869795e+000   -1.436342854982504e+000;...
       -1.048835660976022e+000   -1.115815771583362e+000   -1.184614345408666e+000   -1.256100352867799e+000   -1.331235978799527e+000];
    stblfrac(:,16:20,1) = ...
       [-1.985261982958637e+000   -1.926542865732525e+000   -1.880296841910385e+000   -1.843044812063057e+000   -1.812387604873646e+000;...
       -1.846852935880107e+000   -1.828439745755405e+000   -1.817388844989596e+000   -1.812268962543248e+000   -1.812387604873646e+000;...
       -1.719534615317151e+000   -1.736176665562027e+000   -1.756931455967477e+000   -1.782079727531726e+000   -1.812387604873646e+000;...
       -1.604464355709833e+000   -1.650152416312346e+000   -1.699029550621646e+000   -1.752489822658308e+000   -1.812387604873646e+000;...
       -1.501904088536648e+000   -1.570525854475943e+000   -1.643747672313277e+000   -1.723509779436442e+000   -1.812387604873646e+000;...
       -1.411143947581252e+000   -1.497190629447853e+000   -1.591104422133556e+000   -1.695147748117837e+000   -1.812387604873646e+000];

    stblfrac(:,1:5,2) = ...
       [-4.738866777987500e+002   -1.684460387562537e+001   -5.619926961081743e+000   -3.281734135829228e+000   -2.397479160864619e+000;...
       -2.185953347160669e+001   -3.543320127025984e+000   -1.977029667649595e+000   -1.507632281031653e+000   -1.303310228044346e+000;...
       -2.681009914911080e-001   -4.350930213152404e-001   -5.305212880041126e-001   -6.015232065896753e-001   -6.620641788021128e-001;...
       -9.503065419472154e-002   -1.947070824738389e-001   -2.987136341021804e-001   -3.973064532664002e-001   -4.838698271554803e-001;...
       -1.264483719244014e-001   -2.437377726529247e-001   -3.333750988387906e-001   -4.016893641684894e-001   -4.577316520822721e-001;...
       -1.526287733702501e-001   -2.498255243669921e-001   -3.063859169446500e-001   -3.504924054764082e-001   -3.911254396222550e-001];
    stblfrac(:,6:10,2) = ...
       [-1.959508008521143e+000   -1.708174380583835e+000   -1.550822278332538e+000   -1.447013328833974e+000   -1.376381920471173e+000;...
       -1.199548019673933e+000   -1.144166826374866e+000   -1.115692821970145e+000   -1.103448361903579e+000   -1.101126400280696e+000;...
       -7.174026993828067e-001   -7.694003004766365e-001   -8.178267862332173e-001   -8.615585464741182e-001   -9.003104216523169e-001;...
       -5.579448431371428e-001   -6.215822273361273e-001   -6.771753949313707e-001   -7.267793058476849e-001   -7.720164852674839e-001;...
       -5.069548741156986e-001   -5.523620701546919e-001   -5.956554729327528e-001   -6.378655338388568e-001   -6.796745661620428e-001;...
       -4.309657384679277e-001   -4.709130419301468e-001   -5.113624096299824e-001   -5.525816075847192e-001   -5.948321009341774e-001];
    stblfrac(:,11:15,2) = ...
      [-1.327391983207241e+000   -1.292811209009340e+000   -1.267812588403031e+000   -1.249132310044230e+000   -1.234616432819130e+000;...
       -1.104531584444055e+000   -1.110930462397609e+000   -1.118760810700929e+000   -1.127268239360369e+000   -1.136171639806347e+000;...
       -9.347554970493899e-001   -9.658656088352816e-001   -9.945788535033495e-001   -1.021718797792234e+000   -1.048005562158225e+000;...
       -8.141486817740096e-001   -8.541760575495752e-001   -8.929234555236560e-001   -9.311104141820112e-001   -9.694099704722252e-001;...
       -7.215886443544494e-001   -7.640354693071291e-001   -8.074261467088205e-001   -8.522003643607233e-001   -8.988670244927735e-001;...
       -6.384119892912432e-001   -6.836776839822375e-001   -7.310612144698296e-001   -7.810921001396979e-001   -8.344269070778757e-001];
    stblfrac(:,16:20,2) = ...
       [-1.222879780072203e+000   -1.213041554808853e+000   -1.204541064608597e+000   -1.197016952370690e+000   -1.190232162899989e+000;...
       -1.145449097190615e+000   -1.155224344271089e+000   -1.165719407748303e+000   -1.177246763148178e+000   -1.190232162899989e+000;...
       -1.074094694885961e+000   -1.100624477495892e+000   -1.128270402039747e+000   -1.157812818875688e+000   -1.190232162899989e+000;...
       -1.008502023575024e+000   -1.049129636922346e+000   -1.092166845038550e+000   -1.138712425453996e+000   -1.190232162899989e+000;...
       -9.480479125009214e-001   -1.000533792677121e+000   -1.057363229272293e+000   -1.119941850176443e+000   -1.190232162899989e+000;...
       -8.918931068397437e-001   -9.545526172382969e-001   -1.023797332562095e+000   -1.101496412960141e+000   -1.190232162899989e+000];

    stblfrac(:,1:5,3) = ...
       [-1.354883142615948e+000   -8.855778500552980e-001   -7.773858277863260e-001   -7.357727812399337e-001   -7.181850957003714e-001;...
       -5.193811327974376e-002   -1.633949875159595e-001   -2.617724006156590e-001   -3.392619822712012e-001   -4.018554923458003e-001;...
       -6.335376612981386e-002   -1.297738965263227e-001   -1.985319371835911e-001   -2.624863717000360e-001   -3.174865471926985e-001;...
       -9.460338726038994e-002   -1.756165596280472e-001   -2.282691311262980e-001   -2.638458905915733e-001   -2.918110046315503e-001;...
       -1.158003423724520e-001   -1.620942232133271e-001   -1.790483132028017e-001   -1.937097725890709e-001   -2.109729530977958e-001;...
       -5.695213481951577e-002   -2.485009114767256e-002   -2.455774348005581e-002   -4.243720620421176e-002   -6.906960852184874e-002];
    stblfrac(:,6:10,3) = ...
       [ -7.120493514301658e-001   -7.121454153857569e-001   -7.157018373526386e-001   -7.209253714350538e-001   -7.265425280053609e-001;...
       -4.539746445467862e-001   -4.979328472153985e-001   -5.348184073267474e-001   -5.654705188376931e-001   -5.909430146259388e-001;...
       -3.637544360366539e-001   -4.030045272659678e-001   -4.369896090801292e-001   -4.671253359013797e-001   -4.944847533335236e-001;...
       -3.167744873288179e-001   -3.408290016876749e-001   -3.649204420006245e-001   -3.894754728525021e-001   -4.146904022890949e-001;...
       -2.311198638992638e-001   -2.537077422985343e-001   -2.783252370301364e-001   -3.047045003309861e-001   -3.327092628454751e-001;...
       -1.000745485866474e-001   -1.334091111747126e-001   -1.681287272131953e-001   -2.038409527302062e-001   -2.404547731975402e-001];
    stblfrac(:,11:15,3) = ...
       [-7.317075569303094e-001   -7.359762286696208e-001   -7.392122467978279e-001   -7.414607677550720e-001   -7.428480570989012e-001;...
       -6.123665499489599e-001   -6.307488506465194e-001   -6.469130897780404e-001   -6.615145568123281e-001   -6.750798357120451e-001;...
       -5.198770070249209e-001   -5.439265161390062e-001   -5.671356857543234e-001   -5.899325077218274e-001   -6.127077038151078e-001;...
       -4.406707089221509e-001   -4.675033009839270e-001   -4.952960990683358e-001   -5.242037261193876e-001   -5.544463409264927e-001;...
       -3.623063449447594e-001   -3.935470145089454e-001   -4.265595391976379e-001   -4.615525703717921e-001   -4.988293297210071e-001;...
       -2.780623638274261e-001   -3.168837529800063e-001   -3.572466721186688e-001   -3.995862986780706e-001   -4.444626893956575e-001];
    stblfrac(:,16:20,3) = ...
      [-7.435216571211187e-001   -7.436225251216279e-001   -7.432733099840527e-001   -7.425762029730668e-001   -7.416143171871161e-001;...
       -6.880470899358724e-001   -7.008026232247697e-001   -7.137148222421971e-001   -7.271697520465581e-001   -7.416143171871161e-001;...
       -6.358474023877762e-001   -6.597648782206755e-001   -6.849381555866478e-001   -7.119602076523737e-001   -7.416143171871161e-001;...
       -5.863313160876512e-001   -6.202819599064874e-001   -6.568811178840162e-001   -6.969403639254603e-001   -7.416143171871159e-001;...
       -5.388134824040952e-001   -5.820906647738434e-001   -6.294732446564461e-001   -6.821024214831549e-001   -7.416143171871159e-001;...
       -4.925935308416445e-001   -5.449092276644302e-001   -6.026377433551201e-001   -6.674379829825384e-001   -7.416143171871159e-001];

    stblfrac(:,1:5,4) = ...
       [-4.719005698760254e-003   -5.039419714218448e-002   -1.108600074872916e-001   -1.646393852283324e-001   -2.088895889525075e-001;...
       -3.167687806490741e-002   -6.488347295237770e-002   -9.913854730442322e-002   -1.306663969875579e-001   -1.574578108363950e-001;...
       -6.256908981229170e-002   -1.058190431028687e-001   -1.215669874255146e-001   -1.261149689648148e-001   -1.284283108027729e-001;...
       -7.132464704948761e-002   -5.885471032381771e-002   -3.846810486653290e-002   -2.801768649688129e-002   -2.615407079824540e-002;...
        1.186775035989228e-001    1.847231744541209e-001    1.899666578065291e-001    1.756596652192159e-001    1.538218851318199e-001;...
        1.359937191266603e+000    7.928324704017256e-001    6.068350758065271e-001    4.949176895753282e-001    4.117787224185477e-001];
    stblfrac(:,6:10,4) = ...
       [-2.445873831127209e-001   -2.729819770922066e-001   -2.951510874462016e-001   -3.121233685073350e-001   -3.249196962329062e-001;...
       -1.797875581290475e-001   -1.986122400020671e-001   -2.148458045681510e-001   -2.292024720743768e-001   -2.422125650878785e-001;...
       -1.318108373643454e-001   -1.372885008966837e-001   -1.450218673440198e-001   -1.548461140242879e-001   -1.664940537646226e-001;...
       -3.037902421859952e-002   -3.894619676380785e-002   -5.076849313651704e-002   -6.518223105549245e-002   -8.178056142331483e-002;...
        1.287679439328719e-001    1.022243387982872e-001    7.488543991005173e-002    4.698265181928261e-002    1.852002327642577e-002;...
        3.435869264264112e-001    2.844376471729288e-001    2.312306852681522e-001    1.820841981890349e-001    1.357181057787019e-001];
    stblfrac(:,11:15,4) = ...
       [-3.344714240325961e-001   -3.415532212363377e-001   -3.467713617249639e-001   -3.505859000173167e-001   -3.533413466958321e-001;...
       -2.542699931601989e-001   -2.656748454748664e-001   -2.766656461455947e-001   -2.874428940341864e-001   -2.981872822548070e-001;...
       -1.796994139325742e-001   -1.942454974557965e-001   -2.099854734361004e-001   -2.268483937252861e-001   -2.448403779828917e-001;...
       -1.003134231215546e-001   -1.206343411798188e-001   -1.426762955132322e-001   -1.664453845103147e-001   -1.920257997377931e-001;...
       -1.062008675791458e-002   -4.062891141128176e-002   -7.175196683590498e-002   -1.042870733773311e-001   -1.385948877988075e-001;...
        9.117291945474759e-002    4.766184332000264e-002    4.481886485253039e-003   -3.904933750228177e-002   -8.364689014849616e-002];
    stblfrac(:,16:20,4) = ...
       [-3.552947623689004e-001   -3.566384591258251e-001   -3.575167387322836e-001   -3.580387843935552e-001   -3.582869092425832e-001;...
       -3.090746307371333e-001   -3.202900038682522e-001   -3.320450798333745e-001   -3.445973947956370e-001   -3.582869092425832e-001;...
       -2.640470286750166e-001   -2.846415660837839e-001   -3.069024734642628e-001   -3.312464672828315e-001   -3.582869092425832e-001;...
       -2.195942670864279e-001   -2.494428999135824e-001   -2.820166786810741e-001   -3.179740384308457e-001   -3.582869092425832e-001;...
       -1.751227987938045e-001   -2.144432379167035e-001   -2.573138196343415e-001   -3.047716553689650e-001   -3.582869092425832e-001;...
       -1.301133939768983e-001   -1.794049920724848e-001   -2.327202766583559e-001   -2.916310469293936e-001   -3.582869092425832e-001];

    stblfrac(:,1:5,5) = ...
      [                      0                         0                         0                         0                         0;...
       -2.998229841415443e-002   -3.235136568035350e-002   -1.058934315424071e-002    1.472786013654386e-002    3.649529125352272e-002;...
       -4.911181618214269e-004    7.928758678692660e-002    1.295711243349632e-001    1.575625247967377e-001    1.726794061650541e-001;...
        6.444732609572413e-001    5.412205715497974e-001    4.864603927210872e-001    4.457073928551408e-001    4.118964225372133e-001;...
        4.884639795042095e+000    1.686842470765597e+000    1.132342494635284e+000    8.944978064032267e-001    7.538011200000044e-001;...
        2.410567057697245e+001    4.005534670805399e+000    2.144263118197206e+000    1.518214626927320e+000    1.198109338317733e+000];
    stblfrac(:,6:10,5) = ...
     [                      0                         0                         0                         0                         0;...
        5.320761222262883e-002    6.497369053185199e-002    7.235439352353751e-002    7.603800885095309e-002    7.671459793802817e-002;...
        1.799982238321182e-001    1.821699713013862e-001    1.806145618464317e-001    1.761248753943454e-001    1.691770293512301e-001;...
        3.823074983529713e-001    3.554905959697276e-001    3.305043126978712e-001    3.066571802106021e-001    2.834017043112906e-001;...
        6.558265419066330e-001    5.806408912949470e-001    5.191065509143589e-001    4.663489244354866e-001    4.194539705064985e-001;...
        9.966378800612080e-001    8.532685386168033e-001    7.427048697651345e-001    6.524693172360032e-001    5.756299950589361e-001];
    stblfrac(:,11:15,5) = ...
      [                      0                         0                         0                         0                         0;...
        7.500001602159387e-002    7.139599669434762e-002    6.628276247821394e-002    5.992932695316782e-002    5.250925428603021e-002;...
        1.600901411017374e-001    1.491003610537801e-001    1.363865273697878e-001    1.220722641614886e-001    1.062191001109524e-001;...
        2.602853501366307e-001    2.369238065872132e-001    2.129824521942899e-001    1.881563959610275e-001    1.621474808586950e-001;...
        3.765099860312678e-001    3.361566147323812e-001    2.973499640484341e-001    2.592283952427927e-001    2.210255604589869e-001;...
        5.079606300067100e-001    4.466711396792393e-001    3.897746494263863e-001    3.357416130711989e-001    2.832892169418335e-001];
    stblfrac(:,16:20,5) = ...
      [                      0                         0                         0                         0                         0;...
        4.411421669339249e-002    3.476266163507976e-002    2.439917920106283e-002    1.289010976694223e-002                         0;...
        8.881586460416716e-002    6.976629777350905e-002    4.886974404989612e-002    2.578932638717129e-002                         0;...
        1.346349888095220e-001    1.052403813710735e-001    7.348119932151805e-002    3.870673240105876e-002                         0;...
        1.820030836908522e-001    1.413881485626739e-001    9.829989964989198e-002    5.165115573609639e-002                         0;...
        2.312355801087936e-001    1.783807793433976e-001    1.233869208812706e-001    6.463145748462040e-002    9.714451465470120e-017];


    stblfrac(:,1:5,6) = ...
       [ 4.719005698760275e-003    5.039419714218456e-002    1.108600074872919e-001    1.646393852283322e-001    2.088895889525074e-001;...
        1.944613194060750e-001    3.117984496788369e-001    3.615078716560812e-001    3.879646155737581e-001    4.042606354602197e-001;...
        3.045958300133999e+000    1.315675725057089e+000    9.757973307352019e-001    8.294361410388060e-001    7.456405896421690e-001;...
        2.339312510820383e+001    3.858569195402605e+000    2.091507439545032e+000    1.515362821077606e+000    1.231804842218289e+000;...
        1.231812404655975e+002    9.151933726881032e+000    3.856468345925451e+000    2.470027172456050e+000    1.862167039303084e+000;...
        5.049829135345403e+002    1.890722475322573e+001    6.427275565975617e+000    3.715903402980179e+000    2.636417882085815e+000];
    stblfrac(:,6:10,6) = ...
       [ 2.445873831127209e-001    2.729819770922065e-001    2.951510874462016e-001    3.121233685073347e-001    3.249196962329060e-001;...
        4.152379986226543e-001    4.229018705591941e-001    4.280900470005300e-001    4.311273812611276e-001    4.321442286112657e-001;...
        6.900226415397631e-001    6.495436520935480e-001    6.180526887451320e-001    5.921654464012007e-001    5.697923159645174e-001;...
        1.060749495885882e+000    9.442937075476816e-001    8.583603822642385e-001    7.911221543980916e-001    7.360251815557063e-001;...
        1.521067254392224e+000    1.300039377551776e+000    1.142711537858461e+000    1.023045102736937e+000    9.273664178094935e-001;...
        2.065989355542487e+000    1.711228437455139e+000    1.466088158475343e+000    1.283765226486882e+000    1.140575450959062e+000];
    stblfrac(:,11:15,6) = ...
        [3.344714240325963e-001    3.415532212363379e-001    3.467713617249641e-001    3.505859000173170e-001    3.533413466958320e-001;...
        4.312423594533669e-001    4.285591238013830e-001    4.242644840754073e-001    4.185310514289916e-001    4.115050794489342e-001;...
        5.495326577846258e-001    5.304020801294532e-001    5.116943409858906e-001    4.928954730588648e-001    4.736165965702772e-001;...
        6.890778676134198e-001    6.476526200515113e-001    6.099033923678876e-001    5.744600864566568e-001    5.402514096915735e-001;...
        8.477633920324498e-001    7.792812067953944e-001    7.185943530039393e-001    6.633207377171386e-001    6.116407715135426e-001;...
        1.023262411940948e+000    9.237922892835746e-001    8.369566524681974e-001    7.591595457820644e-001    6.877508180861301e-001];
    stblfrac(:,16:20,6) = ...
       [ 3.552947623689000e-001    3.566384591258254e-001    3.575167387322835e-001    3.580387843935554e-001    3.582869092425831e-001;...
        4.032875933324668e-001    3.939222836649399e-001    3.833860261287606e-001    3.715758694363207e-001    3.582869092425831e-001;...
        4.535361612745278e-001    4.323485980953122e-001    4.097162006469898e-001    3.852184728042033e-001    3.582869092425835e-001;...
        5.063904595668142e-001    4.720865286037160e-001    4.365637761840112e-001    3.989743423180101e-001    3.582869092425835e-001;...
        5.620594176198462e-001    5.132627179036522e-001    4.639774715385669e-001    4.128508865888630e-001    3.582869092425835e-001;...
        6.206265009880273e-001    5.559603894356728e-001    4.919976875425384e-001    4.268552022160075e-001    3.582869092425835e-001];

    stblfrac(:,1:5,7) = ...
       [ 1.354883142615939e+000    8.855778500552969e-001    7.773858277863266e-001    7.357727812399328e-001    7.181850957003700e-001;...
        2.264297017396562e+001    3.703766301758638e+000    2.034998948698223e+000    1.510923485095245e+000    1.265729978744353e+000;...
        1.955956459466261e+002    1.118917023817671e+001    4.357570503031440e+000    2.718083521990130e+000    2.041945502327640e+000;...
        1.131527106972301e+003    2.742019413138009e+001    8.094356141096943e+000    4.405625422851678e+000    3.045873292912599e+000;...
        4.991370610374878e+003    5.832596523112534e+001    1.361736440227531e+001    6.617793943005997e+000    4.277065691957527e+000;...
        1.808482789458792e+004    1.120299053944505e+002    2.131886896428897e+001    9.395528700779570e+000    5.735282952993835e+000];
    stblfrac(:,6:10,7) = ...
       [ 7.120493514301658e-001    7.121454153857567e-001    7.157018373526382e-001    7.209253714350531e-001    7.265425280053608e-001;...
        1.126910935459891e+000    1.039315711942880e+000    9.801156996469297e-001    9.380990288559633e-001    9.070002633955093e-001;...
        1.682687096145072e+000    1.462088170281394e+000    1.313508264506275e+000    1.206803763884095e+000    1.126395471042167e+000;...
        2.368493556832589e+000    1.968378518204384e+000    1.704951233806636e+000    1.518043793772535e+000    1.377948007790416e+000;...
        3.176211386678905e+000    2.549432728119129e+000    2.146593646702069e+000    1.865193645178458e+000    1.656315874739094e+000;...
        4.099439855675913e+000    3.198582996879541e+000    2.632582798272859e+000    2.243339709179312e+000    1.957469852365064e+000];
    stblfrac(:,11:15,7) = ...
       [ 7.317075569303093e-001    7.359762286696208e-001    7.392122467978273e-001    7.414607677550722e-001    7.428480570989009e-001;...
        8.829463516299942e-001    8.633779161543368e-001    8.465599716104961e-001    8.313215935120923e-001    8.168794983145117e-001;...
        1.063360967480519e+000    1.012144436660489e+000    9.690437805764626e-001    9.314651792280744e-001    8.975270882378618e-001;...
        1.268363069256580e+000    1.179563109954373e+000    1.105319244270462e+000    1.041384485194864e+000    9.846979577532636e-001;...
        1.493891969504980e+000    1.362797559741365e+000    1.253624580847262e+000    1.160149469096889e+000    1.078008118654219e+000;...
        1.736744887299007e+000    1.559416515511960e+000    1.412280239489399e+000    1.286729855523644e+000    1.176933895080190e+000];
    stblfrac(:,16:20,7) = ...
       [ 7.435216571211178e-001    7.436225251216276e-001    7.432733099840527e-001    7.425762029730666e-001    7.416143171871158e-001;...
        8.027015701907034e-001    7.884022863227798e-001    7.736657968963813e-001    7.581862145381915e-001    7.416143171871158e-001;...
        8.658237613571567e-001    8.352619776464638e-001    8.049334692839693e-001    7.740056420537431e-001    7.416143171871158e-001;...
        9.329399521299938e-001    8.842632875709708e-001    8.371061471443788e-001    7.900396709438159e-001    7.416143171871157e-001;...
        1.003953952010710e+000    9.354146255148074e-001    8.702022492276336e-001    8.062927602676150e-001    7.416143171871157e-001;...
        1.078670034479511e+000    9.886802003678273e-001    9.042295460529033e-001    8.227686378257326e-001    7.416143171871157e-001];

    stblfrac(:,1:5,8) = ...
        [4.738866777987514e+002    1.684460387562540e+001    5.619926961081758e+000    3.281734135829232e+000    2.397479160864624e+000;...
        4.841681688643794e+003    5.491635522391771e+001    1.256979234254407e+001    6.069209132601843e+000    3.940274296039883e+000;...
        3.154616792561625e+004    1.420805372229245e+002    2.403953052063284e+001    9.998426062380954e+000    5.930362539243756e+000;...
        1.520631636586534e+005    3.148956061770992e+002    4.132943146104890e+001    1.518515134801384e+001    8.367182529059960e+000;...
        5.901656732159231e+005    6.246491282963873e+002    6.581680474603525e+001    2.173557079848703e+001    1.125045444319795e+001;...
        1.944624278667431e+006    1.139848804168331e+003    9.894809619823921e+001    2.974824391888133e+001    1.458002371721213e+001];
    stblfrac(:,6:10,8) = ...
       [1.959508008521145e+000    1.708174380583837e+000    1.550822278332539e+000    1.447013328833976e+000    1.376381920471174e+000;...
        2.963447020215305e+000    2.423693540860402e+000    2.089182215079736e+000    1.865572849084425e+000    1.708118159360888e+000;...
        4.190132768594454e+000    3.268280841745006e+000    2.710662024401290e+000    2.341995909523891e+000    2.082469140437107e+000;...
        5.624308785058203e+000    4.226708866462347e+000    3.402197103627229e+000    2.865360079281767e+000    2.490393899977397e+000;...
        7.254212029229660e+000    5.287806421003054e+000    4.154585933912857e+000    3.428194997160839e+000    2.925780747207696e+000;...
        9.070365685373144e+000    6.442950257298201e+000    4.960971490178073e+000    4.025088868546689e+000    3.384287797654701e+000];
    stblfrac(:,11:15,8) = ...
       [ 1.327391983207241e+000    1.292811209009341e+000    1.267812588403031e+000    1.249132310044230e+000    1.234616432819130e+000;...
        1.593041126030172e+000    1.506471132927683e+000    1.439628954887186e+000    1.386580264484466e+000    1.343153406231364e+000;...
        1.891158929781140e+000    1.745070641877115e+000    1.630251730907927e+000    1.537630629971792e+000    1.460938380853296e+000;...
        2.214464603850502e+000    2.003098342270666e+000    1.835905829230373e+000    1.700021765831942e+000    1.586823477367793e+000;...
        2.557944985263177e+000    2.276562749626175e+000    2.053593165082403e+000    1.871725504345519e+000    1.719630879614922e+000;...
        2.918103805585008e+000    2.562588803694463e+000    2.281050180010934e+000    2.051085944176459e+000    1.858294826115218e+000];
    stblfrac(:,16:20,8) = ...
       [ 1.222879780072204e+000    1.213041554808854e+000    1.204541064608597e+000    1.197016952370690e+000    1.190232162899990e+000;...
        1.306371038922589e+000    1.274091491606534e+000    1.244744203398707e+000    1.217124809801410e+000    1.190232162899990e+000;...
        1.395630981221581e+000    1.338301797693731e+000    1.286320343916442e+000    1.237570697847646e+000    1.190232162899990e+000;...
        1.490188322141933e+000    1.405530485165501e+000    1.329245194088195e+000    1.258353899045780e+000    1.190232162899990e+000;...
        1.589489775546923e+000    1.475587597649461e+000    1.373481210080780e+000    1.279472666002594e+000    1.190232162899990e+000;...
        1.692973560150181e+000    1.548256386823049e+000    1.418980226656540e+000    1.300924242481222e+000    1.190232162899990e+000];

    stblfrac(:,1:5,9) = ...
        [1.890857122067037e+006    1.074884919696010e+003    9.039223076384690e+001    2.645987890965103e+001    1.274134564492299e+001;...
        1.434546473316804e+007    2.987011338973518e+003    1.804473474220022e+002    4.487048929338575e+001    1.960113433547389e+001;...
        7.716266115204613e+007    6.969521346220721e+003    3.196657990381036e+002    6.941784107578008e+001    2.798990029407097e+001;...
        3.253192550565641e+008    1.437315176424486e+004    5.205876769957880e+002    1.006582035946658e+002    3.790739646062081e+001;...
        1.143638705833100e+009    2.703823367877713e+004    7.964291266167923e+002    1.391051003571698e+002    4.935349274736288e+001;...
        3.492208269966229e+009    4.737075925045248e+004    1.161019167208514e+003    1.852377745522907e+002    6.232811767701676e+001];
    stblfrac(:,6:10,9) = ...
        [7.864009406553027e+000    5.591791397752693e+000    4.343949435866960e+000    3.580521076832391e+000    3.077683537175252e+000;...
        1.132727408868559e+001    7.671280872680232e+000    5.732691330034323e+000    4.573075545294608e+000    3.818589092027862e+000;...
        1.533578393202605e+001    9.991349773961725e+000    7.243609507849516e+000    5.634462725204553e+000    4.601857009791827e+000;...
        1.985701175129152e+001    1.252691966593449e+001    8.859346059355138e+000    6.752431092162364e+000    5.418366793527828e+000;...
        2.486500490402286e+001    1.525895955988075e+001    1.056731639206889e+001    7.918478700695184e+000    6.262067266019560e+000;...
        3.033836510475647e+001    1.817240938152932e+001    1.235792736188858e+001    9.126360342186048e+000    7.128676006881803e+000];
    stblfrac(:,11:15,9) = ...
        [2.729262880847459e+000    2.479627528870858e+000    2.297138304998906e+000    2.162196365947915e+000    2.061462692277420e+000;...
        3.297130126832188e+000    2.920640582387343e+000    2.640274592919582e+000    2.426998377788287e+000    2.262233765245289e+000;...
        3.893417077901593e+000    3.382471282597797e+000    2.999860062957988e+000    2.705234908082859e+000    2.473610569743775e+000;...
        4.510891038980249e+000    3.859051363710381e+000    3.370702720510665e+000    2.992693808833481e+000    2.692636527934335e+000;...
        5.144949915764652e+000    4.346635348399592e+000    3.749599176843221e+000    3.286641675099088e+000    2.917178603817272e+000;...
        5.792462636377325e+000    4.842756648977701e+000    4.134472567050430e+000    3.585273662390985e+000    3.145733197974777e+000];
    stblfrac(:,16:20,9) = ...
        [1.985261982958638e+000    1.926542865732524e+000    1.880296841910385e+000    1.843044812063057e+000    1.812387604873647e+000;...
        2.133064562958712e+000    2.029912595114798e+000    1.945516531961286e+000    1.874392545595589e+000    1.812387604873647e+000;...
        2.288441176274372e+000    2.137883347336651e+000    2.012884307837858e+000    1.906295529437326e+000    1.812387604873647e+000;...
        2.449737610939970e+000    2.249772716121334e+000    2.082221357100924e+000    1.938735806854783e+000    1.812387604873647e+000;...
        2.615585030563546e+000    2.364937633815368e+000    2.153342270485199e+000    1.971693892149562e+000    1.812387604873647e+000;...
        2.784907129216124e+000    2.482804054400846e+000    2.226062706102394e+000    2.005149380181030e+000    1.812387604873647e+000];

    %%%% Interpolate to find initial guess
    [alpIn betIn uIn] = meshgrid(alpha,beta,utemp(middle));
    X0(middle) = interp3(Alp,Bet,P,stblfrac,alpIn, betIn, uIn, 'linear');

end

X0 = reshape(X0,size(u));

end





















Contact us