Exploring PDF of Sum of Independent Exponential Random Variables

In this article, we’ll delve into the expressions for the density function of the sum of independent exponential random variables. Additionally, we’ll explore an application involving the probability density function (PDF) of the norm square of a correlated Gaussian random vector.

Table of Contents

Closed-form PDF Expression for the Sum of Independent Exponential Random Variables

Let’s consider $Y$, the sum of $r$ independent exponential random variables, $X_i \sim \exp(\beta_i)$, where $i = 1, 2, \ldots, r$ i.e.,

$$ \begin{align} Y &= X_1 + X_2 + \ldots + X_r \end{align} $$

where $$ \begin{align} f_{X_i}(t) = \beta_i e^{-\beta_i t},\ t \geq 0 \end{align} $$

For simplicity, let’s explore this in three different cases:

Case 1: When all $\beta_i$ values are identical ($\beta_i = \beta$ for all $i$), the PDF of $Y$ is given by:

$$ \begin{align} f_{Y}(t) &= \frac{\beta^r}{(r-1)!}t^{r-1}e^{-t\beta},\ t \geq 0\\ &= {\rm Erl}\left(r,\beta\right)\\ &={\rm Gamma}\left(r,\beta\right),\ r\in \mathbb{Z}_{+} \end{align} $$

This expression resembles Erlang or Gamma distributions for $r \in \mathbb{Z}_{+}$.

Case 2: If all $\beta_i$ values are distinct, the PDF of $Y$ becomes a more intricate expression involving the individual $\beta_i$ values.

\begin{align} f_{Y}(t) &= \sum_{i=1}^r \frac{\beta_1\cdots\beta_r}{\prod_{\substack{j=1\\j\neq i}}\left(\beta_j - \beta_i\right)}e^{-t\beta_i},\ t \geq 0 \end{align}

and

Case 3: In the general scenario, where there are $n$ distinct $\beta$ values with $k_i$ representing the number of times $\beta_i$ is repeated, the density function is provided in the image below:

For detailed derivations, refer to Akkouchi, Mohamed’s work: Akkouchi, Mohamed. “On the convolution of exponential distributions.” Journal of the chungcheong mathematical society 21.4 (2008): 501-501. Link.

Weighted Sum of Exponential Random Variables

Considering a weighted sum of exponential random variables introduces new parameters. $$ \begin{align} Y = \sum_{i=1}^{r}w_i X_i \end{align} $$ The PDF expressions mentioned earlier can be adapted by defining new beta $\bar{\beta}_i = \beta_i/w_i$ where $w_i$ represents the weights of the random variables.

We can also obtain the closed form expression for CDF by making use of lower incomplete gamma function.

Validating Through MATLAB Implementation

You can download my MATLAB implementation of the function that computes PDF/CDF for given input values from MATLAB forumLink. Get the latest version from Github: Link.

The provided script demonstrates the comparison between theoretical PDFs and empirical data histograms for the sum of exponentials. These comparisons serve to validate the derived expressions.

% This script provides an example to demonstrate usage of the function
% MATLAB (tested with 2023a)
% Author: Zakir Hussain Shaik


clc;
clear;
close all;

% Exponential parameters
lmd = [1,1,2,3,4,4];

% Weights of random variables
w = [2,3,4,1,2,5];

% Theoretical expression of PDF
t = 0:0.01:20;
[f,F] = sumexppdf(t,lmd,w);

% Generating emperical samples (N number of samples)
N  = 1e7;

% Weighted sum of exponential random variables
X = 0;
for  i = 1:length(lmd)
    X  = X + w(i)*exprnd(1/lmd(i),N,1);
end


%% Plots
figure;
subplot(2,1,1);
histogram(X,'Normalization','pdf','DisplayStyle','stairs');hold on;
plot(t,f,'*r','MarkerIndices',1:200:length(f));
legend('Emperical','Theoretical','Location','best');
xlabel('$t$','Interpreter','latex');
ylabel('$f_X(t)$','Interpreter','latex');
title('PDF of Sum of Exponentials','Interpreter','latex');
grid on;
xlim([min(t),max(t)]);

subplot(2,1,2);
ecdf(X);hold on;
plot(t,F,'*r','MarkerIndices',1:200:length(f));
legend('Emperical','Theoretical','Location','best');
xlabel('$t$','Interpreter','latex');
ylabel('$F_X(t)$','Interpreter','latex');
title('CDF of Sum of Exponentials','Interpreter','latex');
grid on;
xlim([min(t),max(t)]);

Results/Plots using MATLAB:

This plot showcases the comparison between the theoretical PDF (red markers) and the generated data’s histogram, validating the expressions for the sum of independent Exponentials.

Application: Norm Square of Gaussian Random Vector

Moving beyond the exponential distribution, we explore the squared norm of a complex Gaussian random vector, $\mathbf{z} \sim \mathcal{CN}(\mathbf{0},\mathbf{R})$, where $\mathbf{R} \succ 0$ is a positive definite covariance matrix.

Density of Squared Norm of $\mathbf{z}$ i.e., $\Vert\mathbf{z}\Vert^2$

By leveraging properties of the complex Gaussian vector and employing weighted sums of exponential random variables, we derive the closed-form density function, $f_{\Vert\mathbf{z}\Vert^2}(t)$.

To compute this, let $$ \begin{align} \mathbf{z} = \mathbf{R}^{\frac{1}{2}}\mathbf{x};\ \mathbf{x}\sim \mathcal{CN}(\mathbf{0},\mathbf{I}) \end{align} $$

Then $$ \begin{align} \Vert\mathbf{z}\Vert^2 &= \Vert\mathbf{R}^{\frac{1}{2}}\mathbf{x}\Vert^2\\ &= \Vert\mathbf{U}\boldsymbol{\Sigma}^{\frac{1}{2}}\mathbf{U}^H\mathbf{x}\Vert^2\\ &=\sum_{i=1}^{r}\mu_i \vert y_i\vert^2,\ {y}_i \sim \mathcal{CN}(0,1) \end{align} $$

This is weighted sum of independent central chi-square distributions. Moreover, we know that $\vert y_i\vert^2 \sim \exp(1)$ and therefore we can use weighted sum of exponential random variables where $\beta_i = 1/\mu_i$. Thus, we have the density function, $f_{\Vert\mathbf{z}\Vert^2}(t)$ , in closed form.

This concept finds applications in MIMO systems, where the beamforming gain directly relates to the norm square of the channel vector.

We can also compute PDF of the norm of the Gaussian vector:

Norm of $\mathbf{z}$ i.e, $\Vert\mathbf{z}\Vert$

Furthermore, we compute the PDF of the norm of the Gaussian vector using an expression involving the squared norm.

$$ \begin{align} f_{\Vert\mathbf{z}\Vert}(t) = 2tf_{\Vert\mathbf{z}\Vert^2}(t^2) \end{align} $$

MATLAB Implementations

clc;
clear;
close all;

% Random seed for repeatability
rng(1)

% Dimension of the vector
N = 4;

% Defining a positive definite covariance matrix 
h = crandn(N,N);
R = (h*h' + eye(N));

% Computing eigen values
w = eig(R);

% Theoretical expression of PDF for norm square
t = 0:0.01:120;
lmd = ones(N,1);
f = sumexppdf(t,lmd,w);

% Theoretical expression of PDF for norm
t2 = sqrt(t);
g = 2*sqrt(t).*f;

% random samples
MC = 1e6;
y_normSqr = zeros(MC,1);
y_norm = zeros(MC,1);
for i = 1:MC

    z = sqrtm(R/2)*(randn(N,1) + 1j*randn(N,1));

    y_normSqr(i,1) = norm(z)^2;
    y_norm(i,1) = norm(z);
end



%%
figure;
subplot(2,1,1);
histogram(y_normSqr,'Normalization','pdf','DisplayStyle','stairs','EdgeColor', 'k');hold on;
plot(t,f,'*r','MarkerIndices',1:90:length(f));
legend('Histogram','Theoretical','Location','best');
xlabel('$y$','Interpreter','latex');
ylabel('$f_Y(y)$','Interpreter','latex');
grid on;
title('PDF of $y = \Vert \mathbf{z} \Vert ^2$, $\mathbf{z} \sim \mathcal{CN}(\mathbf{0},\mathbf{R})$','Interpreter','latex');

subplot(2,1,2);
histogram(y_norm,'Normalization','pdf','DisplayStyle','stairs','EdgeColor', 'k');hold on;
plot(t2,g,'*r','MarkerIndices',1:90:length(g));
legend('Histogram','Theoretical','Location','best');
xlabel('$y$','Interpreter','latex');
ylabel('$f_Y(y)$','Interpreter','latex');
grid on;
title('PDF of $y = \Vert \mathbf{z} \Vert$, $\mathbf{z} \sim \mathcal{CN}(\mathbf{0},\mathbf{R})$','Interpreter','latex');

Results/Plots using MATLAB:

Considering Non-Zero Mean

Finally, we delve into a side result involving a Gaussian vector with a non-zero mean. Here, we explore the density of $\Vert\mathbf{z}\Vert^2$ for $\mathbf{z} \sim \mathcal{CN}(\mathbf{m},\mathbf{R})$, where $\mathbf{m}$ represents the non-zero mean.

To compute this, let $$ \begin{align} \mathbf{z} = \mathbf{R}^{\frac{1}{2}}\mathbf{x} + \mathbf{m};\ \mathbf{x}\sim \mathcal{CN}(\mathbf{0},\mathbf{I}) \end{align} $$

Then $$ \begin{align} \Vert\mathbf{z}\Vert^2 &= \Vert\mathbf{R}^{\frac{1}{2}}\left(\mathbf{x} + \mathbf{R}^{-\frac{1}{2}}\mathbf{m}\right)\Vert^2\\ &= \Vert\boldsymbol{\Sigma}^{\frac{1}{2}}\left(\bar{\mathbf{x}} + \mathbf{U}^H\mathbf{R}^{-\frac{1}{2}}\mathbf{m}\right)\Vert^2;\\ &= \Vert\boldsymbol{\Sigma}^{\frac{1}{2}}\left(\bar{\mathbf{x}} + \bar{\mathbf{m}}\right)\Vert^2;\\ &=\sum_{i=1}^{r}\mu_i \vert y_i\vert^2,\ {y}_i \sim \mathcal{CN}(m_i,1) \end{align} $$ where $\bar{\mathbf{x}} = \mathbf{U}^H\mathbf{x}$ (has same distribution as $\mathbf{x}$ because it is circularly symmetric), and $\bar{\mathbf{m}} = \mathbf{U}^H\mathbf{R}^{-\frac{1}{2}}\mathbf{m}$.

This is weighted sum of independent non-central chi-square distributions. We can represent this form as Generalized chi-square distribution:

$$ \begin{align} \Vert\mathbf{z}\Vert^2 &=\sum_{i=1}^{r}\mu_i\vert y_i\vert^2,\ {y}_i \sim \mathcal{CN}(m_i,1)\\ &=\sum_{i=1}^{r}\mu_i\vert\mathcal{R}\left\{y_i\right\}\vert^2 + \vert\mathcal{I}\left\{y_i\right\}\vert^2), \end{align} $$ where

$$ \begin{align} \mathcal{R}\left\{y_i\right\} &\sim \mathcal{N}\left(\mathcal{R}\left\{m_i\right\},\frac{1}{2}\right)\\ \mathcal{I}\left\{y_i\right\} &\sim \mathcal{N}\left(\mathcal{I}\left\{m_i\right\},\frac{1}{2}\right) \end{align} $$

Therefore: $$ \begin{align} \vert\mathcal{R}\left\{y_i\right\}\vert^2 + \vert\mathcal{I}\left\{y_i\right\}\vert^2 &\sim 0.5\chi^2(2,\lambda_i)\\ \lambda_i &= 2\left(\left(\mathcal{R}\left\{m_i\right\}\right)^2 + \left(\mathcal{I}\left\{m_i\right\}\right)^2\right) \end{align} $$

Thus, norm square of Gaussian vector can be modeled as generalized chi-squares as $$ \begin{align} \Vert\mathbf{z}\Vert^2 \sim \widetilde{\chi}^2(\boldsymbol{\omega},\boldsymbol{\kappa},\boldsymbol{\lambda},m,s) \end{align} $$

where $$ \begin{align} \boldsymbol{\omega} &= [\mu_1,\ldots,\mu_r]_{1\times N}\\ \boldsymbol{\kappa}& = [2,\ldots,2]_{1\times N}\\ [\boldsymbol{\lambda}]_i& = 2\left(\left(\mathcal{R}\left\{m_i\right\}\right)^2 + \left(\mathcal{I}\left\{m_i\right\}\right)^2\right)\\ m&=0\\ s & = 0 \end{align} $$ where $[\boldsymbol{\lambda}]_i ,\ i = 1,\ldots,N$ is the $i$th element of vector $[\boldsymbol{\lambda}]$.

Let’s implement this in MATLAB. For this we will be using the helper function from MATLAB File Exchange Central: Link

clc;
clear;

% Random seed for repeatability
rng(1)

% Dimension of the vector
N = 4;

% Defining a positive definite covariance matrix 
%h = sqrt(0.5)*(randn(N,1) + 1j*randn(N,1));
h = crandn(N,N);
R = (h*h' + eye(N));

% Mean
m = sqrt(0.5)*(randn(N,1) + 1j*randn(N,1));

% Emperical samples
MC = 1e6;
y_normSqr = zeros(MC,1);
y_norm = zeros(MC,1);
for i = 1:MC

    z = sqrtm(R/2)*(randn(N,1) + 1j*randn(N,1)) + m;

    y_normSqr(i,1) = norm(z)^2;
    y_norm(i,1) = norm(z);
end

% Computing EVD
[U,S] = eig(R);
S_sqrt = sqrt(S);

% Obtaining modified means
m_b = U'*m;
m_bb = S_sqrt\m_b;

% Set the weights for Generalized Chi-Square
w = diag(S);
k = 2*ones(N,1);

m_r = real(m_bb);
m_i = imag(m_bb);

lambda = m_r.^2 + m_i.^2;

% Theoretical PDF
%Converting all the row dimension
wT = 0.5*w.';
kT = k.';
lambdaT = 2*lambda.';
[f,t]=gx2pdf('full',wT,kT,lambdaT,0,0);

% Collecting only positive values of t
t_pos = t(t>=0);
f_tpos = f(t>=0);

% Theoretical expression of PDF for norm
t2 = sqrt(t_pos);
g = 2*sqrt(t_pos).*f_tpos;

%% Plots
figure;
subplot(2,1,1);
histogram(y_normSqr,'Normalization','pdf','DisplayStyle','stairs','EdgeColor', 'k');hold on;
plot(t_pos,f_tpos,'*r','MarkerIndices',1:600:length(f_tpos));
legend('Histogram','Theoretical','Location','best');
xlabel('$y$','Interpreter','latex');
ylabel('$f_Y(y)$','Interpreter','latex');
grid on;
title('PDF of $y = \Vert \mathbf{z} \Vert ^2$, $\mathbf{z} \sim \mathcal{CN}(\mathbf{m},\mathbf{R})$','Interpreter','latex');

subplot(2,1,2);
histogram(y_norm,'Normalization','pdf','DisplayStyle','stairs','EdgeColor', 'k');hold on;
plot(t2,g,'*r','MarkerIndices',1:600:length(g));
legend('Histogram','Theoretical','Location','best');
xlabel('$y$','Interpreter','latex');
ylabel('$f_Y(y)$','Interpreter','latex');
grid on;
title('PDF of $y = \Vert \mathbf{z} \Vert$, $\mathbf{z} \sim \mathcal{CN}(\mathbf{m},\mathbf{R})$','Interpreter','latex');

Results/Plots using MATLAB:

For real Gaussian random vector, one can use the above methodology to model the norm square with an equivalent generalized-chi-square-distribution with approprtiate parameters and use the function gx2pdf().

For further mathematical insights and detailed derivations related to linear combinations of central/non-central chi-squared distributions:

1. Moschopoulos, P. G., and W. B. Canada. “The distribution function of a linear combination of chi-squares.” Computers & mathematics with applications 10.4-5 (1984): 383-386. Link

2. Mathai, Arakaparampil M., and Serge B. Provost. “Quadratic forms in random variables: theory and applications."(1992). Link