Integral and inverse integral

Question

Mohamed Zied 2023-6-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990148-integral-and-inverse-integral

编辑： Walter Roberson 2024-9-17

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);

3 个评论
显示 1更早的评论隐藏 1更早的评论

Mohamed Zied 2023-6-30

在 MATLAB Online 中打开

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 2023-6-30

syms x

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2023-6-30

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990148-integral-and-inverse-integral#answer_1265473

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Paul 2023-7-3

在 MATLAB Online 中打开

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 2023-7-5

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 2023-7-5

编辑：Paul 2023-7-5

在 MATLAB Online 中打开

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 2023-6-30

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990148-integral-and-inverse-integral#answer_1265383

移动：Torsten 2023-6-30

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Mohamed Zied 2023-7-3

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 2023-7-3

移动：Torsten 2023-7-3

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 2024-9-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990148-integral-and-inverse-integral#answer_1517765

编辑：Walter Roberson 2024-9-17

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2024-9-17

What is the reason that you repeat the question ?

请先登录，再进行评论。

Integral and inverse integral

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Integral and inverse integral

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

2 个评论 显示 无隐藏 无

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

1 个评论
显示 -1更早的评论隐藏 -1更早的评论