Evaluating multivariate stable distributions

3 次查看（过去 30 天）

Mohammad Shojaei Arani 2023-11-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2044747-evaluating-multivariate-stable-distributions

评论： Mohammad Shojaei Arani 2023-11-18

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Mohammad Shojaei Arani 2023-11-10

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 2023-11-10

在 MATLAB Online 中打开

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 2023-11-12

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

Mohammad Shojaei Arani 2023-11-12

Hi Paul,

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

Mohammad Shojaei Arani 2023-11-18

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

请先登录，再进行评论。

请先登录，再回答此问题。