Evaluating multivariate stable distributions

Hello,

I need to calculate the probability density function for the multivariate stable distribution. It seems to me that MATLAB only does so

fo univariate stable distributions (I wish I am wrong). I tested this as below;

pd = makedist('Stable','alpha',0.5,'beta',0.5,'gam',1,'delta',0);

pdf(pd2,1)

ans = 0.1003

pdf(pd2,2)

ans = 0.0533

but when I try the command pdf(pd2,[1 2]) (or pdf(pd2,[1; 2]) ) then I get

0.1003 0.0533

and this means that matlab does not want consider [1 2] as a single point in 2-dimensional state space.

Do you know if there is any package, or anything, where I can perform multivariate estimation of stable densities at some query points?

Thanks for your kind help in advance!

Babak

[Torsten]

The Stable distribution is a one-dimensional distribution. So why should MATLAB consider [1 2] as a point in 2d-space and not as two points in 1d-space ?

[Mohammad Shojaei Arani]

Hi Torsten,

I believe you are wrong. If you take a look at the following wikipedia page you see that it can be multivariate as well:

https://en.wikipedia.org/wiki/Multivariate_stable_distribution

I did not mean that MATLAB should consider [1 2] as a point in 2d (sorry, if I was not clear on this). I meant to say that this was a test to see whether MATLAB interprets [1 2] as point in 2d space or 2 points in 1d space. Unfortunately, the second guess seems to be correct.

[Torsten]

If you look at what MATLAB considers as "Stable Distribution", you'll see that it is univariate:

https://uk.mathworks.com/help/stats/stable-distribution.html

This might help (I didn't look into it more closely):

https://uk.mathworks.com/products/connections/product_detail/stable.html

[Mohammad Shojaei Arani]

Torsten, I am a mathematian and have published papers where I use the concept.

If you say that stable distribution is only univariate then you are wrong (please see the wikipedia page. This means that we must not discuss about this). But, if you say that MATLAB does not calcuate it then it is a different story. However, it is hard to believe that an expensive software like matlab does not have tools to calculate multivariate stable distributions.

[Torsten]

When you type

pd = makedist('Stable','alpha',0.5,'beta',0.5,'gam',1,'delta',0);

you get the univariate Stable distribution.

Unfortunately, I don't see a multivariate extension in MATLAB.

If the Toolbox I gave you the link for does not meet your requirements, I think you will have to program what you want on your own.

What do you want to do with the multivariate stable distribution ?

[Mohammad Shojaei Arani]

Thanks Torsten (sorry I was not polite enough),

Yes. It seems that MATLAB does not have it.

Right now I am concerned with estimating the parameters of stochastic system of differential equations where the stochastic part is Levy noise. Levy noise has stable distribution. Unfortunately, at the moment I can only solve my problem for the case of one-dimensional systems.

But, even for the univariate case I am not going to use the MATLAB built-in functions because it is

computationally very expensive. Fortunately, I could find a MATLAB package by Mike O'Neil which is 100 times faster

(here is the link: https://gitlab.com/s_ament/qastable). Below, gives a comparisson between the 2 methods for the univariate case:

pd = makedist('Stable','alpha',1.6,'beta',0.3,'gam',1,'delta',0);

x=random('Stable',1.6,0.3,1,0,[1 10^5]);

tic;

A=pdf(pd,x); % This is based on MATLAB built in functions

t1=toc;

tic;

B=stable_pdf(x,1.6,0.3); % This is based on the package by Mike O'Neil

t2=toc;

norm(A-B)

t1/t2

ans =

9.3800e-12

ans =

223.9407

[Paul]

@Mohammad Shojaei Arani

Why doesn't stable_pdf require the gamma and delta parameters? Are they assumed to be 1 and 0 respectively? Can they be specified to be other values?

Any idea why stable_pdf is so much faster?

[Mohammad Shojaei Arani]

Hello Torsten,

A univariate stable distribution has 4 parameters (alpha,betta,gamma,delta) where gamma is a scale parameter and somehow plays the role of variance (but it is not vriance since stable distributions might not have a defined variance due to heavy tails). delta is a shift parameter and somehow plays the role of mean. For a standard stable distribution, a gamma =1 and delta=0 will be considered (otherwise, we can always apply a simple linear transformation to normalize it). So, actually the first 2 parameters are very important.

Mike O'Neil has spent a lot of time, in my opinion, to develop optimal codes which are very accurate and fast. Unfortunately, I am not into the details of how he developed such nice codes (my research is different. I just use it as a tool).

But, perhaps you can talk to MATLAB authorities and try to convince them to replace the current MATLAB built-in functions by those of Mike O'Neil. Please note that the package of Mike O'Neil is rather big and can calculate lots of other very useful quntities such as derivatives and integral of the density, etc. In many applications researchers need such quntities and what is important is that his codes are 100's of orders of magnitude faster. Actually, I was going to write a paper. When I realized that it takes 9 seconds to calculate the likelihoods of an array of size 10^5 I was about to give up (note that I need to repeat such a calculation in a for-loop many many times. So, it takes hours if I use built-in functions in MATLAB).

Have a nice day Torsten!

Babak

[Paul]

I think the comment above was to me, not Torsten.

Here is a timing test of Matlab's pdf.

pd = makedist('Stable','alpha',1.6,'beta',0.3,'gam',1,'delta',0);
rng(100);
for count = 1:5
   x = random(pd,[1 10^count]);
   tic;
   A = pdf(pd,x);   % This is based on MATLAB built in functions
   t1(count) = toc;
end
figure
semilogx(10.^(1:5),t1,'-x')

figure
[sx,ind] = sort(x);
plot(sx,A(ind),'-o','MarkerSize',2)
ylim([-0.01, 0.3])

I wonder why the timing jumps up for larger number of values. Maybe it's because for larger number of random values it's more likely that some of the random values will be at points where the pdf is difficult to compute accurately (presumably in the tails?). The matab doc StableDistribution says it computes the pdf via direct integration and perhaps there are some values of x for which it takes more time for the integration to converge. Just speculating ....

Instead of posting norm(A-B), can you post something that compares the actual values of A-B? Maybe a plot of abs(A-B) vs x or abs(A-B)./abs(B) vs x or something like that? I'm curious about how stable_pdf compares pointwise to pdf, particularly in the tails.

[Mohammad Shojaei Arani]

Hello Paul,

Sorry that I am responding rather late (I did not recognize your message before)

Yes. Most probably the timing jumps up for larger values of x. It is not that easy to compte the density of such 'rare events'. I do not know the details but I can intuitively understand this. However, I believe for 'extremely large' values of x

we can benifit from a Theorem (Lemma, whatever) which says that the tails of stable distributions behave as |x|^(-(1+alpha)). But, x should really be big. The following code lines gives the plots you asked me:

pd = makedist('Stable','alpha',1.6,'beta',0.3,'gam',1,'delta',0);

x=random('Stable',1.6,0.3,1,0,[1 10^5]);

tic;

A=pdf(pd,x); % This is based on MATLAB built in functions

t1=toc;

tic;

B=stable_pdf(x,1.6,0.3); % This is based on the package by Mike O'Neil

t2=toc;

norm(A-B)

t1/t2

plot(x,abs(A-B),'-o','MarkerSize',2);

plot(x,abs(A-B)./abs(B),'-o','MarkerSize',2);

I have attached the two figures (I could not figure out how to post (not attach) the figures as you did)

There is just one limitation about the codes of Mike O'Neil I realized recently. It only does the calculations whenever 0.9 < alpha < 1.1 if beta != 0, and alpha < .5 not supported. But, let me tell you that this is never a sad news. In the applications in more than 99% of the cases researchers only need alpha>1, often alpha>1.2 or 1.3. I have never needed to calculate for the cases where alpha<1.4 .

Recently I wrote a code and in my code I use the Mike package whenever his package supports, otherwise I resort to MATLAB. But, my code does not really need MATLAB.

So, I hope MATLAB developers (authorities, whoever) consider this hybrid scheme I told you. In my opinion MATLAB really needs to work on its performance and I often see lots of rooms for the improvement. Some functions, tools, etc are badly slow.

Best,

Babak

[Paul]

For future reference, graphic images, like in a .png file, can be inserted at the cursor by clicking on the image icon (to the left of the link icon in the Insert menu).

Should I be surprised that the largest differences are around x = 0 and not in the tails?

Maybe pdf is slower because it uses the same method to work for all values of alpha and beta (though the doc page does mention an assumption it makes on alpha). Maybe stable_pdf is taking advantage of faster methods that are tuned for its allowable values of parameters, though I suppose there's nothing stopping pdf from doing the same for that same subset of parameters.

If you fell stongly about it, you can always submit an enhancement reqeust to TMW tech support.

openfig('Fig1.fig');

openfig('Fig2.fig');

[Mohammad Shojaei Arani]

Hi Paul,

Now, I learned how to upload figures properly :-)

[Mohammad Shojaei Arani]

Hello Paul,

No. It is not surprising to me that the largest difference occurs arond (not exactly at) x=0 for the following reason reason:

Although the density is heavy tailed but still the density around x=0 is very high. A consequence of this is that when we sample data majority of them are still near 0. If we want to make a fair comparisson we should sample equal amount of data everywhere across the state space. The relative erros abs(A-B)./abs(B) give us a better picture but still it is far from being a 'fair comparisson'.

AUnfortunately, I am extremely busy to investigate this but what I can tell you is that the method of Mike works very well. As I told you, 2 weks ago I started to write a paper and when I noticed that the computational time using MATLAB is very BIG I deciided to give up. Luckilly, I discovered the Mike's codes and write a beautiful paper!

best,

Babak