How to deal with such integration in MATLAB

Hello all, I am working on a problem wherein i have to find analytical expression of outage probability. In my case the outage probability is given as

where are random variables. Specifically, are Gamma random, Z is exponential random and H depends on a factor which is Gaussian random.

Also, the values of constants is as follows: have values such that . Thus if then . ζ can take values like 1, 10, 100, 1000, 10000, 100000 etc. has value as 0.3 and has value like 1.6390,.

I simulated this equation (1) in MATLAB and is working perfectly. But when I tried to obtain analytical expression in terms of X then the problem is that i am getting negative infinity as value of .

The expression that I obtained in terms of X is as follows:

where , and all other terms starting from are PDF of the random variables, is lower incomplete gamma.

I am not getting if equation (1) is working correctly then why equation (2) is giving me negative infinity value.

Any help in this regards will be highly appreciated.

[Torsten]

It would help if you included your MATLAB code in which you plotted a histogram for equation (1).

I guess this is what you mean when you write "I simulated this equation (1) in MATLAB and is working perfectly."

And it would be interesting to know how you theoretically derive (2) from (1).

[chaaru datta]

Thank u sir for ur response. I am sharing the code for simulation of eq. (1) , but note that I cant share a complete code so just sharing the part inside for loop. t1 denotes iterations and we can have it as say 1000.

OP_sim = [];
for zeta = 0:10:40 %zeta is SNR in dB
    zeta_1= 10^(zeta/10); %SNR in linear
    
    
    for t = 1:t1 % number of iterations
        t
        % Generating channels h_0 and h_1
        h_0 = random('Nakagami',m_0,Omega_0); % nakagami rv
        h_0_abs = abs(h_0);
        h_0_abs_sq = (h_0_abs)^2; %gamma rv
        h_1 = random('Nakagami',m_1,Omega_1); % nakagami rv
        h_1_abs = abs(h_1);
        h_1_abs_sq = (h_1_abs)^2; %gamma rv
        Phi_AR = sqrt(variance_Phi_AR_h1 / 2) * (randn + 1j * randn);
        Phi_AR_abs_sq = (abs(Phi_AR))^2;
        h1_eps = sqrt(variance_h1_eps / 2) * (randn + 1j * randn);
        h1_eps_abs_sq = (abs(h1_eps)^2);
        %For user1:
        %for HRIS amplitude response R_1 (eta_1) 
        
        eta_1 = sqrt((sig2_eta_1)/2)*(randn(1,1)+1i*randn(1,1));
        eta_1_abs = abs(eta_1);
        if eta_1_abs <=  Theta_1 
            R1 = A_1; % Assign R_1 = A_1
        else
            R1 = a_11; % Assign R_1 = a_11
        end
        % for HRIS phase response (j_phi_1_tilde)
        phi_1_tilde = 2 * pi * rand; % Generate random phase in the range [0, 2*pi)
        %phi_1_tilde = 0;
        phi_1_phase = exp(1j * phi_1_tilde); %complex exponential e^(j*phi_1_tilde)
        %Final HRIS response (phi_1)
        phi_1 = R1 * phi_1_phase;
        phi_1_abs_sq = (abs(phi_1))^2;
        %For user2:
        %for HRIS amplitude response R_1 (eta_1) 
        eta_2 = sqrt((sig2_eta_2)/2)*(randn(1,1)+1i*randn(1,1));
        eta_2_abs = abs(eta_2);
        if eta_2_abs <=  Theta_2 
            R2 = A_2; % Assign R_2 = A_2
        else
            R2 = a_22; % Assign R_2 = a_2
        end
        % for HRIS phase response (j_phi_2_tilde)
        phi_2_tilde = 2 * pi * rand; % Generate random phase in the range [0, 2*pi)
        %phi_2_tilde = 0;
        phi_2_phase = exp(1j * phi_2_tilde); %complex exponential e^(j*phi_2_tilde)
        %Final HRIS response (phi_2)
        phi_2 = R2 * phi_2_phase;
        phi_2_abs_sq = (abs(phi_2))^2;
        % Parts of outage probability
        Nr = a_1 * zeta_1 * h_0_abs_sq * h_1_abs_sq * g_1_abs_sq * phi_1_abs_sq * P_h1_2;
        Dr_1 = a_2 * zeta_1 * h_0_abs_sq * h_1_abs_sq * g_1_abs_sq * phi_1_abs_sq * P_h1_2;
        Dr_2 = a_1 * zeta_1 * h_0_abs_sq * g_1_abs_sq * phi_1_abs_sq *  Phi_AR_abs_sq;
        
        %Dr_3 = a_1 * zeta_1 * h_0_abs_sq * (P_h1_2*h1_eps_abs_sq) * g_1_abs_sq * phi_1_abs_sq;
        Dr_3 = a_1 * zeta_1 * h_0_abs_sq * (h1_eps_abs_sq) * g_1_abs_sq * phi_1_abs_sq;
        Dr_4 = a_2 * zeta_1 * h_0_abs_sq * g_1_abs_sq* phi_1_abs_sq * Phi_AR_abs_sq;
        
        %Dr_5 = a_2 * zeta_1 * h_0_abs_sq * (P_h1_2*h1_eps_abs_sq) * g_1_abs_sq * phi_1_abs_sq;
        Dr_5 = a_2 * zeta_1 * h_0_abs_sq * (h1_eps_abs_sq) * g_1_abs_sq * phi_1_abs_sq;
        Dr = Dr_1 + Dr_2 + Dr_3 + Dr_4 + Dr_5 + 1;
        
        gamma_1 = Nr/Dr;
        if gamma_1 < u_1
            count = count+1;     
        end  
    end
    OP_sim = [OP_sim,count/(t)];
    count = 0;  
end

Unrecognized function or variable 't1'.

To derive eq.(2) from eq. (1), we follow the following approach:

Rearranging the inequality in eq.(1) in terms of X. Thus we get

Above expression can be written in terms of outage probability as

As X is gamma random variable above expression can be written in terms of CDF of gamma random variable,

To solve above expression we need to average out over and thus we obtain eq. (2).

[Torsten]

I can't believe that you can explicitly derive a CDF for the random variable Nr/Dr (which you seem to claim is the integrand in (2)).

[chaaru datta]

Why sir, did I make any mistake...If so then please let me know.

[Torsten]

The sum of two independent random variables already results in a complicated convolution integral. In the denominator, you have a sum of three different products of random variables - together with the numerator, this cannot result in such a simple threefold integral as written in (2).

[chaaru datta]

I agree with you. But thats why we are averaging with respect to Y, H and Z and so in eq. (2) we are finding three integrals.

[Torsten]

Believe me: your CDF is completely wrong, and you shouldn't make attempts to determine it for such a complicated random variable as given in (1). Each addition and multiplication of random variables will give you a new integral in the expression for the CFD - thus approximately a 9-fold integral would result that can no longer be handled. Stick to your Monte-Carlo method to simulate a large number of samples and use MATLAB's "histogram" with " 'Normalization','pdf' " or " 'Normalization','cdf' " to get an impression of the empirical pdf or cdf.