How to obtain CDF from the below PDF function

I am kind of new to MATLAB and I want to obtain the empitical cumulative distribution function (CDF) of the below PDF:

where:

the PDF function is as follows:

GEV_function=@(y)   (Gamma_d.*Gamma_D)./(Delta_t.*Etta_d.*Etta_D).*(1./y).*(exp(-(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^((-1)./Etta_d))./exp((1-(Mu_D.*Gamma_D)+(Gamma_D.*y)./Delta_t).^((-1)./Etta_D))).*(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^(-1-(1./Etta_d)).*(1-Mu_D.*Gamma_D+Gamma_D.*y./Delta_t).^(-1-(1./Etta_D));
PDF=integral(GEV_function,0,inf);    

To obtain the CDF and to check if I get CDF = 1 as x approuches to infinity, I have integrated the f(x) from minus infinity to positive infinity as follows:

GEV_Original=@(y,x) (Gamma_d.*Gamma_D)./(Delta_t.*Etta_d.*Etta_D).*(1./y).*(exp(-(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^((-1)./Etta_d))./exp((1-(Mu_D.*Gamma_D)+(Gamma_D.*y)./Delta_t).^((-1)./Etta_D))).*(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^(-1-(1./Etta_d)).*(1-Mu_D.*Gamma_D+Gamma_D.*y./Delta_t).^(-1-(1./Etta_D));
CDF = integral2(GEV_Original,0,inf,-inf,inf)

But the answer is 0.6449 and not 1.

I guess the usage of double integral is leading to such an error. Is my code correct?

I appreciate any thoughts or help.

[David Goodmanson]

Hi Manesf,

Mu_d = ? Mu_D = ?

[MANESF]

Oh! Sorry about that. Here are the values for Mu_d and Mu_D:

Moreover here is how and is calculated:

Gamma_D = (gamma(1-Etta_D)-1)/(1-Mu_D)
Gamma_d = (gamma(1-Etta_d)-1)/(1-Mu_d)

I have also updated the original post.

[David Goodmanson]

I ran the integration and got a 'minimum step size reached' error and not .6449. Matlab R2019b update 1.

[MANESF]

I'm using R2016b version whcih gives me 0.6449.

By the way, the values that I provided for the parameters in the original post are just an example. More specifically, , , and are real numbers and can take any value.

[the cyclist]

Can you please double-check what you have posted? I copied your values and formulas directly from here. I then ran the code

x = 0:0.01:1;
y = -1:0.01:1;

[xx,yy] = ndgrid(x,y);


figure
mesh(xx,yy,PDF(xx,yy))

to see what the pdf looked like, and I got negative values:

Here is the code I calculated with. I don't think I mistranscribed anything.

Delta_t = 25;
Etta_d = 0.02;
Etta_D = 0.08;
Mu_d = -2.95;
Mu_D = 0.05;
Gamma_d = (gamma(1-Etta_d) - 1)./(1-Mu_d);
Gamma_D = (gamma(1-Etta_D) - 1)./(1-Mu_D);


PDF=@(x,y) (Gamma_d.*Gamma_D)/(Delta_t.*Etta_d.*Etta_D).*(1./y).*(exp(-(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^((-1)./Etta_d))./exp((1-(Mu_D.*Gamma_D)+(Gamma_D.*y)./Delta_t).^((-1)./Etta_D))).*(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^(-1-(1./Etta_d)).*(1-Mu_D.*Gamma_D+Gamma_D.*y./Delta_t).^(-1-(1./Etta_D));

[Jesus Sanchez]

Isn't this a problem of the order of the integral operations? First the PDF should be integrated in terms of y and then the result should be integrated in terms of x to get the CDF

[MANESF]

@the cyclist, thanks a lot for your reply. My bad. I made a mistake in reporting the PDF in the original code as an integral is missing in my original post. Sorry about that. Here is the correct code:

GEV_function=@(y) (1./y).*(exp(-(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^((-1)./Etta_d))./exp((1-(Mu_D.*Gamma_D)+(Gamma_D.*y)./Delta_t).^((-1)./Etta_D))).*(1-(Mu_d.*Gamma_d)+(Gamma_d.*x)./y).^(-1-(1./Etta_d)).*(1-Mu_D.*Gamma_D+Gamma_D.*y./Delta_t).^(-1-(1./Etta_D));
PDF=(Gamma_d*Gamma_D)/(Delta_t*Etta_d*Etta_D)*integral(GEV_function,0,inf);

I have also editted the original post to reflect this change.

The code is now consistant with the formula in the original post as follows:

[MANESF]

@Jesus Sanchez, thanks for your reply. That's why I used double intrgral to deal with it. It might be wrong though, I'm not sure, but I am not able to do it in MATLAB in two seperate integrals (1-integrating PDF with regards to y and 2-integrating the results with regard to x) as the first integral also contains x and MATLAB requires the integrant to be a numerical value. So it is not possible to have a symbolic form of the PDF in terms of x without an integral part for the second integration.

[David Goodmanson]

Hi Manesf,

So far it looks like x ( -inf<=x<=inf ) and y ( 0 <= y <= inf ) can be varied completely independently. In that case, suppose that 1/Etta_d is not an integer. There is sure to be a region of x and y where

( 1 - Mu_d.*Gamma_d  + Gamma_d *x/y )

is negative, in which case

( 1 - Mu_d.*Gamma_d  + Gamma_d *x/y ) ^( (-1)/Etta_d) )

is complex, not real. So it appears that something is not quite right.

Answer 1

[David Goodmanson]

Hello Manesf,

This pdf seems related to the Weibull distribution. I made some abrreviations

Etta_d = xid       Etta_D = xiD
Gamma_d = gd       Gamma_D = gD
Mu_d = mud         Mu_D = muD,

and changed the sign of the argument of the second exp, multiplying by it rather than dividing. This leads to

PDF = @(x,y) (gd*gD)/(Del*xid*xiD)*(1./y) ...
 .*exp(-(1 -mud*gd +gd*x./y ).^(-1/xid)).*(1 -mud*gd +gd*x./y ).^(-1-1/xid) ...   [1]
 .*exp(-(1 -muD*gD +gD*y/Del).^(-1/xiD)).*(1 -muD*gD +gD*y/Del).^(-1-1/xiD);      [2]

Take a look at line [2] , considered as a distribution on its own. There is a basic normalized distribution (whose name I don't know),

(1/xiD) * exp(-y.^(-1/xiD)) .* y.^(-1-1/xiD)

with 0 <= y <= inf. Variable y cannot be negative since it is taken to a fractional power.

Line [2] contains a scaled and shifted version of y. Now the argument itself must be positive,

0 <= (1 -muD*gD +gD*y/Del) <= inf 

The net result is that y can take on some negative values, and the mean value of y turns out to be y_mean = 1, which I assume is intentional. The lower limit for y is

y0 = ((muD*gD-1)*Del/gD)

which is negative. If you integrate line [2] from y0 to inf, and include the normalization factors in front that are associated with that line,

gD/(Del*xiD)  

you get 1.

For the full pdf, integration over y is complcated because now you have the condition in line [1]

0 <= (1 -mud*gd +gd*x./y ) <= inf

and the y limits depend on x. But you sensibly wanted to do both the x integration and the y integration as a check, which is fairly easy. Substitute

x = y*z dx = y*dz

and do the z integration first. The y in y*dz cancels the 1/y in the pdf and the resulting integration is from z0 to inf, where z0 can be calculated in the same way that y0 was. The result is

Del = 25;
xid  = .02;
xiD = .08;
mud = -2.95;
muD = .05;
gD = (gamma(1-xiD)-1)/(1-muD);
gd = (gamma(1-xid)-1)/(1-mud); 
PDF1 = @(z,y) (gd*gD)/(Del*xid*xiD) ...
 .*exp(-(1 -mud*gd +gd*z    ).^(-1/xid)).*(1 -mud*gd +gd*z    ).^(-1-1/xid) ...
 .*exp(-(1 -muD*gD +gD*y/Del).^(-1/xiD)).*(1 -muD*gD +gD*y/Del).^(-1-1/xiD);
z0 = ((mud*gd-1)/gd)*.9999;
y0 = ((muD*gD-1)*Del/gD)*.9999;
CDF1 = integral2(PDF1,z0,inf,y0,inf)
format long
CDF1 =
   0.999999898129661

  

There are some integration issues near z0 and y0. To address that, I just multiplied each limit by .9999. To really do it right would take some work but I think the simple approach here shows that the result is correct.

[MANESF]

Hi David, Thanks a lot for taking the time to answer my question in detail. Much appreciated!