MATLAB Answers

How to obtain CDF from the below PDF function

40 views (last 30 days)
MANESF
MANESF on 16 Dec 2019
Commented: MANESF on 18 Dec 2019
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.

  9 Comments

Show 6 older comments
MANESF
MANESF on 17 Dec 2019
@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
MANESF on 17 Dec 2019
@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
David Goodmanson on 17 Dec 2019
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.

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 18 Dec 2019
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.

  1 Comment

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

Sign in to comment.

More Answers (0)

Sign in to answer this question.


Translated by