Trying to calculate a function with several integrals

Hi,

I'm trying to calculate the difference between two random variables which are following different disribution laws. Then, I try to calculate the probability that the difference is superior to 0. Here is the code :

syms x y z f(x) g(y) a(y,z) c(z)

sigmaR = 10 ; sigmaC = 25 ; muR = 450 ; muC = 390 ; % paramètres des lois normales
k = 40 ; lambda = 479.6; teta = 0; %paramètre de la loi de Weibull
AR = 1/(sigmaR*sqrt(2*pi()));
AC = 1/(sigmaC*sqrt(2*pi()));


f(x) = AC .* exp((-1/2).*((x-muC)./sigmaC).^2) ; % Loi normale 1
g(y) = AR .* exp((-1/2).*((y-muR)./sigmaR).^2) ; % Loi normale 2
%g(y) = (k/lambda) * (y/lambda)^(k-1) * exp(-(y/lambda)^k); % Loi de Weibull


a(y,z) = f(z+y)*g(y);
c(z) = int(a,y,[-Inf +Inf]); %Densité de probabilité de X-Y

double(c(1))
P = (int(c,z,[0 +Inf])) %Probabilité que C soit > 0

hold on
grid on
fplot(f,[200 600])
fplot(g,[200 600])
fplot(c,[-30 30])
% u=[-200:1:30];
% d=u;
% for i=1:length(u)
%     d(i)=c(u(i));
% end
% plot(u,d);
hold off

When i'm executing the code with g(y) being the normal law 2, it works well. But, when I try with g(y) = weibull law, it doesn't work :

• the function c(z) computes well because the code line "double(C(1))" gives a coherent value,
• The P variable doesn't gives a value but a function,
• The curve fplot(c,[-30 30]) is completely false (it does not match the result of plot(u,d) which is in commentary

I've tried to transform the code to use "function_handle" instead of symbolic but without success.

I hope that someone will find how to resolve this issue. Thank you in advance.

Answer 1

[Torsten]

You won't get an analytic expression for the density of the convolution of the normal distribution and the Weilbull distribution.

I suggest to use Monte-Carlo simulation to determine the probability for X-Y>0.

Another option would be numerical integration.

sigmaR = 10 ; sigmaC = 25 ; muR = 450 ; muC = 390 ; % paramètres des lois normales
k = 40 ; lambda = 479.6; teta = 0; %paramètre de la loi de Weibull
n = 1e8;
X = muC + randn(n,1)*sigmaC; 
Y1 = muR + randn(n,1)*sigmaR;
Y2 = wblrnd(lambda,k,n,1);
npos1 = nnz(X-Y1>0)

npos1 = 1292500

npos1/n

ans = 0.0129

npos2 = nnz(X-Y2>0)

npos2 = 378489

npos2/n

ans = 0.0038

[Pierre]

Hello Torsten,

Thanks a lot for your answer, it works perfectly using Monte-Carlo. You also suggest to use numerical integration, I'm interested that you develop this suggestion because I tried something like this but without success and maybe the numerical integration can help me when probabilities are very low because using Monte-Carlo doesn't work when I want to calculate probabilities littler than 1/n :

sigmaR = 10 ; sigmaC = 25 ; muR = 450 ; muC = 390 ; % paramètres des lois normales
AR = 1/(sigmaR*sqrt(2*pi));
AC = 1/(sigmaC*sqrt(2*pi));
k = 40 ; lambda = 479.6; teta = 0; %paramètre de la loi de Weibull
f = @(x) AC .* exp((-1/2).*((x-muC)./sigmaC).^2) ;
g = @(y) AR .* exp((-1/2).*((y-muR)./sigmaR).^2) ;
%g = @(y) (k/lambda) * (y/lambda)^(k-1) * exp(-(y/lambda)^k);

a = @(y,z) f(z+y).*g(y) ;
c = @(z) integral(@(y) a (z,y),-Inf,+Inf) ;
P = integral(c,0,+Inf) 

Error using integralCalc/finalInputChecks
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Thanks again.

[Torsten]

The results don't look promising. I seems difficult for the integrator to correctly reproduce where - given z - f(z+y)*g(y) differs significantly from 0.

sigmaR = 10 ; sigmaC = 25 ; muR = 450 ; muC = 390 ; % paramètres des lois normales
AR = 1/(sigmaR*sqrt(2*pi));
AC = 1/(sigmaC*sqrt(2*pi));
k = 40 ; lambda = 479.6; teta = 0; %paramètre de la loi de Weibull
f = @(x) AC * exp((-1/2).*((x-muC)./sigmaC).^2) ;
g = @(y) AR * exp((-1/2).*((y-muR)./sigmaR).^2) ;
%g = @(y) (k/lambda) * (y/lambda).^(k-1) .* exp(-(y/lambda).^k);
c = @(z)integral(@(y)f(z+y).*g(y),-Inf,Inf,'ArrayValued',true);
z = -1000:0.1:1000;
plot(z,c(z))

% This should be 1, I guess
integral(@(z)c(z),-Inf,Inf,'ArrayValued',true)

ans = 1.9872e-41

P = integral(@(z)integral(@(y)f(z+y).*g(y),-Inf,Inf,'ArrayValued',true,'RelTol',1e-14,'AbsTol',1e-14),0,Inf,'ArrayValued',true,'RelTol',1e-14,'AbsTol',1e-14)

P = 5.9666e-53

[Pierre]

Oups, it seems that this problem is really difficult to solve! Thank you Torsten for spending your time to answer this question. Even if the solution is not the one I expected because of the limitation of Monte-Carlo with very low probabilities, it is still a solution!

I don't know if I let this question open, maybe someone will find another way to make the numerical integration works.

Thank's again Torsten.