c语言 exponent alpha,beta;答案,SaS稳定分布函数（symmetric D -stable distribution）

function r = stblrnd(alpha,beta,gamma,delta,varargin)

%STBLRND alpha-stable random number generator.

% R = STBLRND(ALPHA,BETA,GAMMA,DELTA) draws a sample from the Levy

% alpha-stable distribution with characteristic exponent ALPHA,

% skewness BETA, scale parameter GAMMA and location parameter DELTA.

% ALPHA,BETA,GAMMA and DELTA must be scalars which fall in the following

% ranges :

% 0 < ALPHA <= 2

% -1 <= BETA <= 1

% 0 < GAMMA < inf

% -inf < DELTA < inf

%

%

% R = STBLRND(ALPHA,BETA,GAMMA,DELTA,M,N,...) or

% R = STBLRND(ALPHA,BETA,GAMMA,DELTA,[M,N,...]) returns an M-by-N-by-...

% array.

%

%

% References:

% [1] J.M. Chambers, C.L. Mallows and B.W. Stuck (1976)

% "A Method for Simulating Stable Random Variables"

% JASA, Vol. 71, No. 354. pages 340-344

%

% [2] Aleksander Weron and Rafal Weron (1995)

% "Computer Simulation of Levy alpha-Stable Variables and Processes"

% Lec. Notes in Physics, 457, pages 379-392

%

if nargin < 4

error('stats:stblrnd:TooFewInputs','Requires at least four input arguments.');

end

% Check parameters

if alpha <= 0 || alpha > 2 || ~isscalar(alpha)

error('stats:stblrnd:BadInputs',' "alpha" must be a scalar which lies in the interval (0,2]');

end

if abs(beta) > 1 || ~isscalar(beta)

error('stats:stblrnd:BadInputs',' "beta" must be a scalar which lies in the interval [-1,1]');

end

if gamma < 0 || ~isscalar(gamma)

error('stats:stblrnd:BadInputs',' "gamma" must be a non-negative scalar');

end

if ~isscalar(delta)

error('stats:stblrnd:BadInputs',' "delta" must be a scalar');

end

% Get output size

[err, sizeOut] = genOutsize(4,alpha,beta,gamma,delta,varargin{:});

if err > 0

error('stats:stblrnd:InputSizeMismatch','Size information is inconsistent.');

end

%---Generate sample----

% See if parameters reduce to a special case, if so be quick, if not

% perform general algorithm

if alpha == 2 % Gaussian distribution

r = sqrt(2) * randn(sizeOut);

elseif alpha==1 && beta == 0 % Cauchy distribution

r = tan( pi/2 * (2*rand(sizeOut) - 1) );

elseif alpha == .5 && abs(beta) == 1 % Levy distribution (a.k.a. Pearson V)

r = beta ./ randn(sizeOut).^2;

elseif beta == 0 % Symmetric alpha-stable

V = pi/2 * (2*rand(sizeOut) - 1);

W = -log(rand(sizeOut));

r = sin(alpha * V) ./ ( cos(V).^(1/alpha) ) .* ...

( cos( V.*(1-alpha) ) ./ W ).^( (1-alpha)/alpha );

elseif alpha ~= 1 % General case, alpha not 1

V = pi/2 * (2*rand(sizeOut) - 1);

W = - log( rand(sizeOut) );

const = beta * tan(pi*alpha/2);

B = atan( const );

S = (1 + const * const).^(1/(2*alpha));

r = S * sin( alpha*V + B ) ./ ( cos(V) ).^(1/alpha) .* ...

( cos( (1-alpha) * V - B ) ./ W ).^((1-alpha)/alpha);

else % General case, alpha = 1

V = pi/2 * (2*rand(sizeOut) - 1);

W = - log( rand(sizeOut) );

piover2 = pi/2;

sclshftV = piover2 + beta * V ;

r = 1/piover2 * ( sclshftV .* tan(V) - ...

beta * log( (piover2 * W .* cos(V) ) ./ sclshftV ) );

end

% Scale and shift

if alpha ~= 1

r = gamma * r + delta;

else

r = gamma * r + (2/pi) * beta * gamma * log(gamma) + delta;

end

end

%==== function to find output size ======%

function [err, commonSize, numElements] = genOutsize(nparams,varargin)

try

tmp = 0;

for argnum = 1:nparams

tmp = tmp + varargin{argnum};

end

if nargin > nparams+1

tmp = tmp + zeros(varargin{nparams+1:end});

end

err = 0;

commonSize = size(tmp);

numElements = numel(tmp);

catch

err = 1;

commonSize = [];

numElements = 0;

end

end