How to obtain CDF from the below PDF function

조회 수: 28 (최근 30일)
MANESF
MANESF 2019년 12월 16일
댓글: MANESF 2019년 12월 18일
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
MANESF
MANESF 2019년 12월 17일
@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 2019년 12월 17일
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
David Goodmanson 2019년 12월 18일
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
MANESF
MANESF 2019년 12월 18일
Hi David, Thanks a lot for taking the time to answer my question in detail. Much appreciated!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by