Doubt in solving nested integral

I have a query regarding MATLAB code to plot SOP using analytical expression given by :

where

\begin{equation}

\begin{aligned}

\Theta(u,y,v)

&= \frac{1}{A_k(u,y,v)}

\exp\!\left(

\frac{\lambda_Y \lambda_Z}{\alpha A_k(u,y,v)}

- A_k(u,y,v)\frac{1}{\bar{\gamma}} \right. \\

&\quad \left.

- \lambda_U u - \lambda_V v

\right)

E_1\!\left(

\frac{\lambda_Y \lambda_Z}{\alpha A_k(u,y,v)}

\right).

\end{aligned}

\end{equation}

Here is what i tried for code of SOP (analytical):

% This code is for analytical expression of SOP of user k
clc; clear all; close all;
%% Parameters
lambda_X = 1; % rate parameter of the exponential distribution
lambda_Y = 1;
lambda_Z = 1;
lambda_U = 1;
lambda_V = 1;
alpha = 0.1;  
k = 1;
a1 = 0.1;
a2= 0.3; 
a3 = 0.6;
a = [a1 a2 a3];
R1 = 0.1;   % nearest
R2 = 0.2;   % middle
R3 = 0.3;   % farthest
R = [R1 R2 R3];
beta = 0.8;
t1 = 1*1e5; % samples
t2 = 1e4; 
D = (1 - beta) + beta * sum(a(1:k-1));
SOP_ana = [];
for snr_user_dB = 0:5:80
    
    snr_user = 10^(snr_user_dB/10);
    snr_eve_dB = 0; % in dB
    snr_eve = 10^(snr_eve_dB/10);
    accum = 0;
    for iter = 1:t1
        
        % Generate all RVs together
        y = exprnd(1/lambda_Y);
        u = exprnd(1/lambda_U);
        v = exprnd(1/lambda_V);
        % Gamma_e
        Gamma_e = (beta * a(k) * u) / ...
            (u*D + alpha*y*v + (1/snr_eve));
        % t(u,y,v)
        t_val = (2^R(k))*(1 + Gamma_e) - 1;
        % A_k
        if beta*a(k) - t_val*D <= 0
            continue;
        end
        A_k = (lambda_X * t_val) / (beta*a(k) - t_val*D);
        % Theta
        arg = (lambda_Y * lambda_Z)/(alpha*A_k);
        Theta = (1/A_k) * ...
            exp(arg - A_k*(1/snr_user) ...
            - lambda_U*u - lambda_V*v) ...
            * expint(arg);
        accum = accum + Theta;
    end
    SOP = 1 - (lambda_U * lambda_V * lambda_Y * lambda_Z / alpha) ...
        * (accum / t1);
    SOP_ana = [SOP_ana SOP];
end
snr_axis = 0:5:80;
figure;
%semilogy(snr_axis, SOP_ana, 'r-o','LineWidth',1.5);
plot(snr_axis, SOP_ana, 'r-o','LineWidth',1.5);
grid on;
xlabel('SNR (dB)');
ylabel('SOP');
title('SOP vs SNR');

The output of above code should match the simulation of SOP whose code is given by

% This code is for simulation expression of SOP of user k
%clc; clear all; close all;
%% Parameters
lambda_X = 1; % rate parameter of the exponential distribution
lambda_Y = 1;
lambda_Z = 1;
lambda_U = 1;
lambda_V = 1;
alpha = 0.1;  
k = 1;
a1 = 0.1;
a2= 0.3; 
a3 = 0.6;
a = [a1 a2 a3];
R1 = 0.1;   % nearest
R2 = 0.2;   % middle
R3 = 0.3;   % farthest
R = [R1 R2 R3];
beta = 0.8;
t1 = 10*1e4; % samples
D = (1 - beta) + beta * sum(a(1:k-1));
SOP_sim = []; count = 0;
for snr_user_dB = 0:5:80
    
    snr_user = 10^(snr_user_dB/10);
    snr_eve_dB = 0; % in dB
    snr_eve = 10^(snr_eve_dB/10);
    
    
    for iter = 1:t1
        X = exprnd(1/lambda_X); % mean is 1/lambda
        Y = exprnd(1/lambda_Y);
        Z = exprnd(1/lambda_Z);
        U = exprnd(1/lambda_U);
        V = exprnd(1/lambda_V);
        
        Gamma_k = (beta * a(k)  * X) / ...
            ( X * D + alpha * Y * Z + (1/snr_user));
        Gamma_e =  (beta * a(k)  * U) / ...
            ( U * D + alpha * Y * V + (1/snr_eve));
        
        if Gamma_k < ((2^R(k))*(1+Gamma_e)-1)
            %if Gamma_k < (2^R(k)-1)
            count = count + 1;
        end 
        
        
    end
    
    SOP = count / t1;
    
    SOP_sim = [SOP_sim SOP];
    count = 0;
end
snr_axis = 0:5:80;
figure;
%semilogy(snr_axis, SOP_sim, 'r-o','LineWidth',1.5);
plot(snr_axis, SOP_sim, 'r-o','LineWidth',1.5);
grid on;
xlabel('SNR (dB)');
ylabel('SOP');
title('SOP vs SNR');

I am not getting what wrong i am doing in analytical SOP code...

Any help in this regard is highly appreciated

[Torsten]

@charu shree

Is this interpretation of the missing LaTeX part correct ?

Otherwise, please change it using the LaTeX button in the INSERT section of the menu bar.

I don't understand your Monte-Carlo integration over [0,Inf]^3. Usually, you take a box [0,M]^3, sample N uniformly distributed triples (x1i,x2i,x3i) from this box and evaluate M^3 * 1/N * sum_{i=1}^{i=N} f(x1i,x2i,x3i) to approximate your three-fold nested integral. Could you explain your method ?

[Walter Roberson]

@Torsten

I did not know how to handle the \left and \right or the \quad when I was converting the \begin{equation} form, so I left that block as raw latex instead of converting it.