Re: Levy ('Stable') Distributions in Octave

[Cédric Févotte]

Hi Steven,

The following piece of code generates samples from the alpha-stable 
distribution, of which the Levy distribution is a particular case 
(alpha=1/2,beta=1), with pdf given there

http://en.wikipedia.org/wiki/L%C3%A9vy_distribution

The code generates values for c=1 (scale parameter).

It requires sampling from a Gamma distribution
p(x) \propto x^(alpha-1) * exp(-beta*x)

and then uses the Chambers Mallows Stuck method.

If you don't have a function for generating gamma samples let me know.

Best regards



function X=starnd(alpha,beta,N)

% STARND draws from standard stable distributions
% (Chambers-Mallow-Stuck method).
%
% X=starnd(alpha,beta,N)
%
% Inputs:
% - alpha is the characteristic exponent (0<alpha<=2)
% - beta is skewness parameter (-1<beta<1)
% - N is the number of samples
%
% Output:
% - X is a 1 x N vector containing the sampled values
%
% Linked m-files: this function uses the gamrnd.m Statistics Matlab Toolbox
% function for sampling from Gamma distribution.
%
% Author: C. FEVOTTE
% Email: cedricfevotte AT free DOT fr
% (please report any bug)
% Reference: R. Weron, "On the Chambers-Mallows-Stuck method for
% simulating skewed stable random variables", Technical Report, Hugo
% Steinhaus Center for Stochastic Methods, Wroclaw, 1996

if nargin==2;
    N=1;
end

V=rand(1,N)*pi-pi/2;
W=gamrnd(ones(1,N),ones(1,N));

if alpha==1
    X=(2/pi)*((pi/2+beta*V).*tan(V)-beta*log(W.*cos(V)./(pi/2+beta*V)));
else
   B=atan(beta*tan(pi*alpha/2))/alpha;
   S=(1+beta^2*tan(pi*alpha/2)^2)^(1/(2*alpha));

X=S.*sin(alpha*(V+B))./(cos(V).^(1/alpha)).*(cos(V-alpha*(V+B))./W).^((1-alpha)/alpha);
end

Pav, Steven E wrote:
> Hi;
>
> I cannot find an implementation of Levy Distributions (pdf, cdf, inv,
> rnd) in octave or octave-forge, or on the internet (google search of
> 'stable distribution octave' comes up with all kinds of links for the
> other meanings of stable and distribution.) Did I just not see it?
> It seems like sampling could be just a wrapper for the GSL version, but
> the pdf and cdf might be tricky.  There is a toolbox for this in Matlab
> (200 bucks, www.robustanalysis.com).  Any hints?
>
>
> thanks,
>  
>
> --sep
>
>
> [ Steven E. Pav, Ph.D.            address@hidden ]
> [ Senior Scientist, Tech. Dev.               WORK: 925.463.4253 ]
> [ Nellcor, Tyco Healthcare                   FACS: 925.463.4210 ]
> [ a palindrome:                               tracker trek cart ]
> [ Please Note: The information in this E-mail message may be    ]
> [ legally privileged and confidential information intended only ]
> [ for the use of the individual(s) named above. If you, the     ]
> [ reader of this message, are not the intended recipient, you   ]
> [ are hereby notified that you should not further disseminate,  ]
> [ distribute, or forward this E-mail message. If you have       ]
> [ received this E-mail in error, please notify the sender and   ]
> [ promptly delete the original message. Thank you.              ]
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www.cae.wisc.edu/mailman/listinfo/help-octave

--
Dr Cédric Févotte
Researcher
Mist Technologies
204 rue de Crimée, 75019 Paris, France
Tel: +33 1 55 26 42 31 - Fax: +33 1 55 26 42 03
Web: http://persos.mist-technologies.com/~cfevotte