# How to obtain CDF from the below PDF function

40 views (last 30 days)
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.

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 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 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.

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) ... 
.*exp(-(1 -muD*gD +gD*y/Del).^(-1/xiD)).*(1 -muD*gD +gD*y/Del).^(-1-1/xiD); 
Take a look at line  , 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  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  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 
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 on 18 Dec 2019
Hi David, Thanks a lot for taking the time to answer my question in detail. Much appreciated!