Coding Random Variable Generation in MATLAB, Python & C++

Aug 30, 2023 13 min read Simulations

Welcome to this blog post, where we will learn to generate random variable (RV) using three popular programming languages: MATLAB, Python, and C++. Our focus will be on various types of random variables, including uniform, Gaussian, and Laplacian distributions, or any other distribution with a known inverse Cumulative Distribution Function (CDF). Don’t worry about optimization – the code examples provided prioritize clarity. This tutorial presumes readers have basic coding background, atleast in MATLAB.

Let’s kick things off by examining the Probability Density Functions (PDFs) of these random variables:

Uniform distribution from $a$ to $b$ $f (x) = {\begin{cases} \frac{1}{b - a}, & if a \leq x \leq b \\ 0, & otherwise \end{cases}$

Gaussian distribution $f (x) = \frac{1}{σ \sqrt{2 π}} \cdot e^{- \frac{(x - μ)^{2}}{2 σ^{2}}}$

Laplace distribution $f (x) = \frac{1}{2 b} \cdot e^{- \frac{| x - μ |}{b}}$ where $μ$ is a location parameter and $b > 0$ is a scale parameter.

Let’s begins with MATLAB. For uniform and Gaussian random variables, MATLAB offers rand and randn functions, respectively. However, as of 2023, there’s no built-in function for generating Laplace random variables. Fear not, though; the MATLAB Exchange community provides code solutions. In this post, we’ll explore the inverse-sampling theorem to generate Laplace distribution samples. For that we need inverse CDF of Laplace distribution which is given as following

$F^{- 1} (u) = μ - b sign (u - 0.5) \log (1 - 2 | u - 0.5 |)$

In order to generate samples from a given RV, such as $Y$ , for which the inverse Cumulative Distribution Function $F^{- 1} (y)$ is known, we setps according to inverse-sampling theorem. First, a sample $u$ is drawn from a uniform distribution over the interval $[0, 1)$ . This $u$ is then used to derive a sample $y$ from the RV $Y$ , where $y = F^{- 1} (u)$ . This two-step approach generate samples from various distributions with known inverse CDFs.

In Python, the numpy library’s random module comes to the rescue. You can employ numpy.random.rand() for uniform RVs, numpy.random.randn() for Gaussian RVs, and numpy.random.laplace() (or the inverse-sampling theorem approach) for Laplace RVs.

For C++, you can harness the power of the default_random_engine class template for random number generation and the uniform_real_distribution and normal_distribution classes for uniform and Gaussian distributions, respectively. While no dedicated class exists for Laplace distribution up to C++11, we can utilize the inverse-sampling theorem method. Alternatively, the IT++ library streamlines the process with its range of functionalities (IT++ library have Laplace RVs generating function). IT++ is a C++ library for linear algebra, signal processing, communications, and other related mathematical operations. It provides a wide range of classes and functions for various tasks like matrix operations, signal generation, modulation, error correction coding, and more.

Necessary code explanations are provided after each code section. While our examples will involve generating RVs with specified mean and variance, remember that certain distributions, like the Poisson or Exponential, mean and variance are inter dependent.

For given mean $μ$ and variance $σ^{2}$ , the parameters of uniform distribution are $a = μ - \sqrt{3} σ$ and $b = μ + \sqrt{3} σ$ , and for Laplace distribution $b = \frac{σ}{\sqrt{2}}$ .

Table of Contents

Generate RVs using MATLAB

%{
 This Matlab script generates : Uniform, Gaussian and Laplacian RVs
 License: This code is licensed under the GPLv2 license.

 Author: Zakir Hussain Shaik
%}

clc;
clearvars;
close all;

% Number of samples to generate for each RV
N = 1e6;

% Mean and Variance of RV
Xmean = 0;
Xvar  = 1;

% Generate Uniform RV

% Giving lower and upper limits
a = Xmean - sqrt(3*Xvar);
b = Xmean + sqrt(3*Xvar);

% Generate N samples of uniform distribution
X1 = a + (b-a).*rand(N,1);

% Theoretical PDF of the generated X1 RV
% range for visualization -r:r
r = 5; % number r is just for visual purpose
n1 = a-r:0.01:b+r; 
X1_th = (1/(b-a)).*(n1>a & n1<b);

% Generate Gaussian RV
X2 = Xmean + sqrt(Xvar).*randn(N,1);

% Theoretical PDF of the generated X1 RV
n2 = min(X2):0.01:max(X2);
X2_th = (1/sqrt(2*pi*Xvar))*exp(-((n2 - Xmean).^2)/(2*Xvar));

% Generate Laplacian RV

% Inverse sampling theorem
% First obtain the inverse CDF of Laplacian distribution
% Then plug-in uniform RV (0,1) into that function
u = rand(N,1);
X3 = Xmean -sqrt(Xvar/2)*sign(u-0.5).*log(1-2*abs(u-0.5));

% Theoretical PDF of the generated X3 RV
n3 = min(X3):0.01:max(X3);
X3_th = (1./(sqrt(2*Xvar))).*exp(-abs(n3 - Xmean)/sqrt(Xvar/2));

%% All the required plots
% Plotting Uniform RV
figure;
histogram(X1,'Normalization','pdf'); hold on;
plot(n1,X1_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
dim = [0.2 0.5 0.3 0.3];
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Uniform random variable');

% Plotting Gaussian RV
figure;
histogram(X2,'Normalization','pdf'); hold on;
plot(n2,X2_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Gaussian random variable');

% Plotting Laplacian RV
figure;
histogram(X3,'Normalization','pdf'); hold on;
plot(n3,X3_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Laplacian random variable');

In the above MATLAB code: rand and randn generates a sample from uniform distribution from $[0, 1)$ and normal Gaussian random variable, respectively. One can generate vector samples as done in the code.

Generate RVs using Python

"""
This python script generates : Uniform, Gaussian and Laplacian RVs
License: This code is licensed under the GPLv2 license.

@author: Zakir Hussain Shaik
"""

# Import numpy and matplotlib libraries
import numpy as np
import matplotlib.pyplot as plt

# Number of samples to generate for each RV
N = int(1e6)

# Mean of RV
Xmean = 0
Xvar  = 1

# Generate Uniform RV
a = Xmean - np.sqrt(3*Xvar)
b = Xmean + np.sqrt(3*Xvar)
X1 = a + (b-a)*np.random.rand(N)

# range for visualization -r:r
r = 5 # number r is just for visual purpose
n1 = np.arange(a-r, b+r, 0.01) 
X1_th = (1/(b-a))*((n1>a) & (n1<b))

# Generate Gaussian RV
X2 = Xmean + np.sqrt(Xvar)*np.random.randn(N)

n2 = np.arange(min(X2), max(X2), 0.01)
X2_th = (1/np.sqrt(2*np.pi*Xvar))*np.exp(-((n2 - Xmean)**2)/(2*Xvar))

# Generate Laplacian RV

# Inverse sampling theorem
# First obtain the inverse CDF of Laplacian distribution
# Then plug-in uniform RV (0,1) into that function
# u = np.random.rand(N)
# Xmean -np.sqrt(Xvar/2)*np.sign(u-0.5)*np.log(1-2*np.abs(u-0.5))

# Directly using python function
X3 =np.random.laplace(Xmean,np.sqrt(Xvar/2),N)

n3 = np.arange(min(X3), max(X3), 0.01)
X3_th = (1./(np.sqrt(2*Xvar)))*np.exp(-np.abs(n3 - Xmean)/np.sqrt(Xvar/2))

## All the required plots
# Plotting Uniform RV
plt.figure()
plt.hist(X1, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n1, X1_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Uniform random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Plotting Gaussian RV
plt.figure()
plt.hist(X2, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n2, X2_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Gaussian random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Plotting Laplacian RV
plt.figure()
plt.hist(X3, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n3, X3_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Laplacian random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Show all the plots
plt.show()

In the above Python code: np.random.rand(), np.random.randn() and np.random.laplace() generates a sample from uniform distribution from $[a, b)$ , Gaussian random variable with mean $μ$ and variance $σ^{2}$ , and Laplace distribution with parameters $μ$ and $b$ respectively. One can generate vector samples as done in the code.

Generate RVs using C++11

#define _USE_MATH_DEFINES

#include <iostream>
using namespace std;

#include <fstream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <random>
int main() {
    // Number of samples to generate for each RV
    int N = 1e6;

    // Mean of RV
    double Xmean = 0;
    double Xvar = 1;

    // Generate Uniform RV
    double a = Xmean - sqrt(3 * Xvar);
    double b = Xmean + sqrt(3 * Xvar);
    vector<double> X1(N);
    default_random_engine genU_RV1;
    uniform_real_distribution<double> dist1(a, b);
    
    for (int i = 0; i < N; ++i) {
        X1[i] = dist1(genU_RV1);
    }

    double r = 5;
    vector<double> n1;
    for (double val = a - r; val <= b + r; val += 0.01) {
        n1.push_back(val);
    }
    vector<double> X1_th(n1.size());
    for (size_t i = 0; i < n1.size(); ++i) {
        X1_th[i] = (1 / (b - a)) * (n1[i] > a && n1[i] < b);
    }

    // Generate Gaussian RV
    default_random_engine genG_RV;
    normal_distribution<double> distribution(Xmean, sqrt(Xvar));
    vector<double> X2(N);
    for (int i = 0; i < N; ++i) {
        X2[i] = distribution(genG_RV);
    }

    vector<double> n2;
    for (double val = *min_element(X2.begin(), X2.end()); val <= *max_element(X2.begin(), X2.end()); val += 0.01) {
        n2.push_back(val);
    }
    vector<double> X2_th(n2.size());
    for (size_t i = 0; i < n2.size(); ++i) {
        X2_th[i] = (1 / sqrt(2 * M_PI * Xvar)) * exp(-pow(n2[i] - Xmean, 2) / (2 * Xvar));
    }

    // Generate Laplacian RV
    default_random_engine genU_RV2;
    uniform_real_distribution<double> dist2(0, 1);
    vector<double> X3(N);
    for (int i = 0; i < N; ++i) {
        double u = dist2(genU_RV2);
        X3[i] = Xmean - sqrt(Xvar / 2) * (u < 0.5 ? -log(2 * u) : log(2 - 2 * u));
    }

    vector<double> n3;
    for (double val = *min_element(X3.begin(), X3.end()); val <= *max_element(X3.begin(), X3.end()); val += 0.01) {
        n3.push_back(val);
    }
    vector<double> X3_th(n3.size());
    for (size_t i = 0; i < n3.size(); ++i) {
        X3_th[i] = (1 / (sqrt(2 * Xvar))) * exp(-(abs(n3[i] - Xmean) / sqrt(Xvar / 2))); // careful with abs. abs does absolute of integer. So it should precede with std:: (if not used using namespace)  to overide with function in cmath library
    }
    
    // Save to CSV files
    ofstream file1("X1.csv");
    for (size_t i = 0; i < X1.size(); ++i) {
        file1 << X1[i] << "\n";
    }
    file1.close();

    ofstream file2("X2.csv");
    for (size_t i = 0; i < X2.size(); ++i) {
        file2 << X2[i] << "\n";
    }
    file2.close();

    ofstream file3("X3.csv");
    for (size_t i = 0; i < X3.size(); ++i) {
        file3 << X3[i] << "\n";
    }
    file3.close();

    ofstream file_n1_X1_th("n1_X1_th.csv");
    for (size_t i = 0; i < n1.size(); ++i) {
        file_n1_X1_th << n1[i] << "," << X1_th[i] << "\n";
    }
    file_n1_X1_th.close();

    ofstream file_n2_X2_th("n2_X2_th.csv");
    for (size_t i = 0; i < n2.size(); ++i) {
        file_n2_X2_th << n2[i] << "," << X2_th[i] << "\n";
    }
    file_n2_X2_th.close();

    ofstream file_n3_X3_th("n3_X3_th.csv");
    for (size_t i = 0; i < n3.size(); ++i) {
        file_n3_X3_th << n3[i] << "," << X3_th[i] << "\n";
    }
    file_n3_X3_th.close();

    return 0;
}

In the above C++ code: default_random_engine is a random number engine class from the C++ standard library’s <random> header that generates pseudo-random numbers, genU_RV1 is an instance of the default_random_engine. Next uniform_real_distribution is another class template from the header. It represents a distribution of real numbers that are uniformly distributed within a specified range. dist1 is an instance of the uniform_real_distribution class and is used to define the distribution properties. The line X1 = dist1(genU_RV1); generates a random double-precision floating-point number using the random number generator engine genU_RV1 and the distribution dist1, and assigns the generated value to the variable X1. Similar explanation holds for other distributions. These results are then saved in .csv file to analyze/plot them in either MATLAB or using Python.

Generate RVs using C++11 using IT++ library

#define _USE_MATH_DEFINES

#include <iostream>
#include <fstream>
#include <itpp/itbase.h>
#include <itpp/stat/misc_stat.h>

int main() {
     // Initialize random number generator
    itpp::RNG_randomize();

    // Number of samples to generate for each RV
    int N = 1e6;

    // Mean of RV
    double Xmean = 0;
    double Xvar = 1;

    // Generate Uniform RV
    itpp::Uniform_RNG uniform_rng(Xmean - sqrt(3 * Xvar), Xmean + sqrt(3 * Xvar));
    itpp::vec X1 = uniform_rng.operator()(N);

    // Generate Gaussian RV
    itpp::vec X2 = itpp::randn(N);

    // Generate Laplacian RV
    itpp::Laplace_RNG laplace_rng(Xmean, Xvar);
    itpp::vec X3 = laplace_rng.operator()(N);

    // Save to CSV files
    double a = Xmean - sqrt(3 * Xvar);
    double b = Xmean + sqrt(3 * Xvar);
    double r = 5.0; // Threshold range for visualization
    itpp::vec n1 = itpp::linspace(a-r,b+r,1000);
    itpp::vec X1_th(n1.size());
    for(int i=0;i<n1.size();++i)
    {
      if(n1[i]>a && n1[i]<b)
      {
        X1_th[i] = 1.0 / (b-a);
      }
      else
      {
        X1_th[i] = 0;
      }
    }

    itpp::vec n2 = itpp::linspace(itpp::min(X2), itpp::max(X2), 1000);
    itpp::vec X2_th = itpp::exp(-itpp::pow(n2 - Xmean,2) / (2 * Xvar)) / sqrt(2 * M_PI * Xvar);

    itpp::vec n3 = itpp::linspace(itpp::min(X3), itpp::max(X3), 1000);
    itpp::vec X3_th = itpp::exp(-itpp::abs(n3 - Xmean) / sqrt(Xvar / 2)) / (sqrt(2 * Xvar));

    //std::cout <<itpp::min(X3)<< ", "<<itpp::max(X3) <<", n3[500]=" <<n3[500] << ", X3_th[500]="<< X3_th[500]<< std::endl;
	// Save to CSV files
    std::ofstream file1("X1.csv");
    for (size_t i = 0; i < X1.size(); ++i) {
        file1 << X1[i] << "\n";
    }
    file1.close();

    std::ofstream file2("X2.csv");
    for (size_t i = 0; i < X2.size(); ++i) {
        file2 << X2[i] << "\n";
    }
    file2.close();

    std::ofstream file3("X3.csv");
    for (size_t i = 0; i < X3.size(); ++i) {
        file3 << X3[i] << "\n";
    }
    file3.close();

    std::ofstream file_n1_X1_th("n1_X1_th.csv");
    for (size_t i = 0; i < n1.size(); ++i) {
        file_n1_X1_th << n1[i] << "," << X1_th[i] << "\n";
    }
    file_n1_X1_th.close();

    std::ofstream file_n2_X2_th("n2_X2_th.csv");
    for (size_t i = 0; i < n2.size(); ++i) {
        file_n2_X2_th << n2[i] << "," << X2_th[i] << "\n";
    }
    file_n2_X2_th.close();

    std::ofstream file_n3_X3_th("n3_X3_th.csv");
    for (size_t i = 0; i < n3.size(); ++i) {
        file_n3_X3_th << n3[i] << "," << X3_th[i] << "\n";
    }
    file_n3_X3_th.close();


    return 0;
}

In the above C++ code with IT++ library: Uniform_RNG is a uniform random variable constructor from the IT++ standard library whose inputs are used to set $a$ and $b$ values i.e, limits of uniform random variable. In uniform_rng.operator() , operator() is a member function to draw corresponding distribution random samples. For Gaussian one can use constructor approach or directly use the function randn from IT++ library. For Laplace random variable samples generation, the explanation similar to uniform random variable holds, mutatis mutandis. The function linspace is same as that in MATLAB.

Plotting Curves in MATLAB Using Data from .csv Files Generated by C++ Code

%{
 This Matlab script generates : Uniform, Gaussian and Laplacian RVs
 License: This code is licensed under the GPLv2 license.

 Author: Zakir Hussain Shaik
%} 

clc;
clearvars;
close all;

% Mean of RV
Xmean = 0;
Xvar  = 1;

% Reading .cpp code generated data from .csv
P1 = readmatrix('n1_X1_th.csv'); n1 = P1(:,1); X1_th = P1(:,2);
P2 = readmatrix('n2_X2_th.csv'); n2 = P2(:,1); X2_th = P2(:,2);
P3 = readmatrix('n3_X3_th.csv'); n3 = P3(:,1); X3_th = P3(:,2);

% Reading Random variables
X1 = readmatrix('X1.csv');
X2 = readmatrix('X2.csv');
X3 = readmatrix('X3.csv');

%% All the required plots
% Plotting Uniform RV
figure;
histogram(X1,'Normalization','pdf'); hold on;
plot(n1,X1_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
dim = [0.2 0.5 0.3 0.3];
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Uniform random variable');

% Plotting Gaussian RV
figure;
histogram(X2,'Normalization','pdf'); hold on;
plot(n2,X2_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Gaussian random variable');

% Plotting Laplacian RV
figure;
histogram(X3,'Normalization','pdf'); hold on;
plot(n3,X3_th,'r','LineWidth',2);
legend('Histogram','Theoretical','Location','best');
str = {['Mean = ',num2str(Xmean)],['Variance = ',num2str(Xvar)]};
annotation('textbox',dim,'String',str,'FitBoxToText','on');
xlabel('x');
ylabel('f_X(x)');
title('PDF: Laplacian random variable');

Plotting Curves in Python Using Data from .csv Files Generated by C++ Code

"""
This python script generates : Uniform, Gaussian and Laplacian RVs
License: This code is licensed under the GPLv2 license.

@author: Zakir Hussain Shaik
"""

import numpy as np
import matplotlib.pyplot as plt

# Mean and variance of random variable
Xmean = 0
Xvar = 1

# Reading .csv generated data
P1 = np.loadtxt('n1_X1_th.csv', delimiter=',')
n1, X1_th = P1[:, 0], P1[:, 1]

P2 = np.loadtxt('n2_X2_th.csv', delimiter=',')
n2, X2_th = P2[:, 0], P2[:, 1]

P3 = np.loadtxt('n3_X3_th.csv', delimiter=',')
n3, X3_th = P3[:, 0], P3[:, 1]

# Reading random variables
X1 = np.loadtxt('X1.csv', delimiter=',')
X2 = np.loadtxt('X2.csv', delimiter=',')
X3 = np.loadtxt('X3.csv', delimiter=',')

plt.figure()
plt.hist(X1, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n1, X1_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Uniform random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Plotting Gaussian RV
plt.figure()
plt.hist(X2, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n2, X2_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Gaussian random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Plotting Laplacian RV
plt.figure()
plt.hist(X3, density=True, bins=200) # Histogram with normalized frequency
plt.plot(n3, X3_th, 'r', linewidth=2) # Theoretical PDF
plt.legend(['Theoretical', 'Histogram'], loc='best')
plt.xlabel('x')
plt.ylabel('f_X(x)')
plt.title('PDF: Laplacian random variable')
plt.text(0.2, 0.5, f'Mean = {Xmean}\nVariance = {Xvar}', transform=plt.gca().transAxes) # Text box with mean and variance

# Show all the plots
plt.show()