Integral and inverse integral

Hi,

I would like to calculate the probability of failure using the convolution therom as described in this paper https://www.researchgate.net/publication/314278481_Reliability_Index_for_Non-normal_Distribution_of_Limit_State_Functions

First I wanted to write the code using examples of this paper (Please look at the attached screenshots). Unfortunately, I didn't get the same results. I found a probabilty of failure Pf=3.82 10^-08 instead of 5.55 10^-02.

In addition, I am struggling to write the inverse function to derive the reliabilty index. I always leads to errors?

Can anyone explain the mistake?

Sigma= 1;
Mu=5; 
PDF_Norm=@(x) exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*pi));
a=2;
b=1;
Gamma=@(x) x.^(a-1).*exp(-x/b)/((b^a).*gamma(a));
FUN=@(x) PDF_Norm(x).*Gamma(-x); 
Pf=integral(@(x) FUN(x),-Inf,0);
sym x
Beta=finverse(F,x);

[Torsten]

Where do you define F from the command

Beta=finverse(F,x);

?

[Mohamed Zied]

Actually, it is just a typing error.

Here is the right one. But even that way, x is not recognized.

Unrecognized function or variable 'x'.

Error in untitled (line 10)

Beta=finverse(FUN,x);

sym x
Beta=finverse(FUN,x);

[Walter Roberson]

syms x

Answer 1

[Paul]

The attachment convolution.jpg is confusing. Equation (15) defines Pf(x) as a convolution. But the unnumbered equation in the right column defines Pf as number that is the integral (over x, presumbably) of a convolution. The expression for beta at the bottom of the right column makes no sense.

syms x
Sigma= 1;
Mu=5; 
PDF_Norm(x) = exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*pi));
a=2;
b=1;
Gamma(x) = x.^(a-1).*exp(-x/b)/((b^a).*gamma(a));

% compute Pf(x) (a function of x) as the convolution R(x)*Q(-x) (equation 15 in attachment)
syms w
Pf(x) = simplify(int(Gamma(-w)*PDF_Norm(x-w),w,-inf,inf))

Pf(x) = 

% compute Pf (a number) as the integral of R(x)*Q(-x)
Pf = int(Pf(x),x,-inf,0)

Pf = 

vpa(Pf) % matches attachment

ans = 

0.05554498269121153825992156462484

[Mohamed Zied]

Thank you for your answer, it helps a lot.

You said that the reliability index expression (beta) makes no sense (mathematically, I guess). Could you explain more about it? By definition, the reliability index is the negative inverse of the cumulative density function of the limit state function, Pf in my case. How can I overcome this issue if Matlab is not able to find a functional inverse? can I find the numerical value without finding the functional inverse?

Thank you in advance.

Beta=finverse(Pf,x);

Warning: Unable to find functional inverse.

[Paul]

In the expressions for beta, both equation (16) and the example on the lower right side, the integrand is the convolution R(x)*Q(-x), which is a function of x. But the variable of integration is tau, which makes no sense (at least to me).

In equation (16) I also don't understand the R>Q notation, but I'm sure that's just because I'm not familiar with the material.

If beta is the functional inverse of another function, then beta itself should be a function, but the right side of the page shows that beta is a number (similar to how Pf is treated both as a function and as a number).

The code becomes a lot cleaner by using symbolic constants throughout. And, I multiplied Gamma(x) by heaviside because the gamma pdf is 0 for x < 0.

syms x
Sigma = sym(1);
Mu    = sym(5); 
PDF_Norm(x) = exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*sym(pi)));
a = sym(2);
b = sym(1);
Gamma(x) = x.^(a-1).*exp(-x/b)/((b^a).*gamma(a))*heaviside(x);
% compute Pf(x) (a function of x) as the convolution R(x)*Q(-x) (equation 15 in attachment)
syms w
Pf(x) = simplify(int(Gamma(-w)*PDF_Norm(x-w),w,-inf,inf))

Pf(x) = 

Pfval = vpa(int(Pf(x),x,-inf,0))

Pfval = 

0.055544996883718716151702333141112

Here's the plot of Pf(x), which is a probability density function of the random variable X.

figure
fplot(Pf(x),[-5 10])

And here the cumulative distribution function of X (since the cumulative density function was mentioned in the comment).

Cf(x) = int(Pf(w),w,-inf,x)

Cf(x) = 

figure
fplot(Cf(x),[-5 10]),grid

But I'm still not sure how, or even if, Cf(x) is supposed to relate to beta.

[Mohamed Zied]

By definition, the reliability index is the negative inverse of the CDF of the limit state function. In most cases, the limit state function is G=R-Q; where R is the resistance (Normal distribution in our case) and Q is the load (follows Gamma distribution in our case).

[Paul]

Hi Mohamed,

I think you meant to say: By definition, the reliability index is the negative of the inverse of the CDF of the limit state function evaluated at the probability of failure.

That definition, which seems to be what the authors are saying in the text above equation (13) and in equation (13) itself, is incorrect.

Let G = R - Q with F_G(x) defined as the CDF of G, i.e., F_G(x) = P(G < x).

By defintion, the probability of failure, P_F, is:

P_F = P(G < 0) (equation (3)) = F_G(0) (by definition of the CDF of G).

Therefore, it must be true, always, that 0 = invF_G(P_F).

We can easily show this using the code from above

syms x real
Sigma = sym(1);
Mu    = sym(5); 
PDF_Norm(x) = exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*sym(pi)));
a = sym(2);
b = sym(1);
Gamma(x) = x.^(a-1).*exp(-x/b)/((b^a).*gamma(a))*heaviside(x);

% compute the pdf of G (equation 15)
syms w
P_f(x) = simplify(int(Gamma(-w)*PDF_Norm(x-w),w,-inf,inf));

% compute the probability of failure
P_F = int(P_f(x),x,-inf,0);

% compute the CDF of G
F_G(x) = int(P_f(w),w,-inf,x);

figure
fplot(F_G(x),[-5 10]),grid
xlabel('x'); ylabel('F_G(x) = P(G < x)')

% graphically, the value of the inverse of F_G(x) evaluated at P_F is found
% by drawing a horizontal line at P_F and picking the value off the x-axis
% at the point of intersection with F_G(x), which is at x = 0, as must be
% the case based on how P_F is defined
yline(double(P_F))

Also, the stated definition doesn't work in the case when R and Q are both normal. In this case, the reliability index, beta, is given by equation (1), where Phi(x) is the inverse of the CDF of the standard normal distribution, which easily follows from the definition of P_F. When R and Q are normal, G, the limit state function, is normal, but it is not, in general, standard normal, e.g., as shown in example 1.

Unforturnately, the authors just state equation (16) (which is notationally suspect) without giving any context as to how that equation was derived or what it's supposed to mean (as is shown in Figure 3 for the normal-normal case). It kind of looks like it's supposed to be an implementation of equation (13). But, as discussed above, equation (13), as written, looks incorrect, and I'm not sure that (16) actually implements (13) anyway. So I'm not sure what's going on with equation (16).

Disclaimer: I'm not a structural engineer.

Answer 2

[Torsten]

MATLAB is not able to find the inverse:

syms x
Sigma= 1;
Mu=5; 
PDF_Norm=exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*pi));
a=2;
b=1;
Gamma= x.^(a-1).*exp(-x/b)/((b^a).*gamma(a));
FUN= PDF_Norm*subs(Gamma,x,-x)

FUN = 

Pf=int(FUN,x,-Inf,0)   

Pf = 

Beta=finverse(FUN)

Warning: Unable to find functional inverse.

 
Beta =
 
Empty sym: 0-by-1
 

[Mohamed Zied]

Thank you for your answer, it is helpful.

But is there a way to derive the negative inverse of Pf (as shown in the attached screenshots)?

Thank you in advance.

[Torsten]

Given y, the inverse function of Pf is defined by the equation Pf(x) - y = 0.

To numerically determine x, use a root finder, e.g. "fzero" or "fsolve".

Answer 3

[Sonali]

syms x
Sigma= 1;
Mu=5; 
PDF_Norm=exp(-0.5.*((x-Mu)/Sigma).^2)/(Sigma*sqrt(2*pi));
a=2;
b=1;
Gamma= x.^(a-1).*exp(-x/b)/((b^a).*gamma(a));
FUN= PDF_Norm*subs(Gamma,x,-x)
Pf=int(FUN,x,-Inf,0)   
Beta=finverse(FUN)

Warning: Unable to find functional inverse.

 
Beta =
 
Empty sym: 0-by-1
 

Warning: Unable to find functional inverse.

Beta =

Empty sym: 0-by-1

[Torsten]

What is the reason that you repeat the question ?