Python, MATLAB, Julia, R code: Chapter 4

Acknowledgement: The Julia code is written by the contributors listed here.

Chapter 4.3

CDF of a Gaussian mixture

% MATLAB code to generate the PDF and CDF
gmm = gmdistribution([0; 5], cat(3,1,1), [0.3 0.7]);

x  = linspace(-5, 10, 1000)';
f  = pdf(gmm, x);
x1 = linspace(-5,1,1000)';
f1 = pdf(gmm,x1);
figure;
area(x1,f1,'EdgeColor','none','FaceColor',[0.8,0.8,1]);hold on;
plot(x, f, 'LineWidth', 6);
grid on;
axis([-5 10 0 0.4]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);

F = cdf(gmm, x);
figure;
plot(x, F, 'LineWidth', 6);
grid on;
axis([-5 10 0 1]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);
# Julia code to generate the PDF and CDF
using Distributions
using Plots

gmm = MixtureModel(Normal[Normal(0, 1), Normal(5, 1),], [0.3, 0.7])

p1 = plot(x -> 0, -5, 1, fillrange= x->pdf(gmm, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false) # style
plot!(p1, x -> pdf(gmm, x), -5, 10,
    color=1, linewidth=6, label="PDF") # style

p2 = plot(x -> cdf(gmm, x), -5, 10,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

CDF of a uniform random variable

% MATLAB code to generate the PDF and CDF of a uniform random variable
unif = makedist('Uniform','lower',-3,'upper',4);

x  = linspace(-5, 10, 1500)';
f  = pdf(unif, x);
x1 = linspace(-5,1,1500)';
f1 = pdf(unif,x1);
figure;
area(x1,f1,'EdgeColor','none','FaceColor',[0.8,0.8,1]);hold on;
plot(x, f, 'LineWidth', 6);
grid on;
axis([-5 10 0 0.4]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);

F = cdf(unif, x);
figure;
plot(x, F, 'LineWidth', 6); hold on;
plot(x(600), F(600), 'ko', 'LineWidth', 2, 'MarkerSize', 10, 'MarkerFaceColor', [0.8,0.8,1]);
grid on;
axis([-5 10 0 1]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);
# Python code to generate the PDF and CDF of a uniform random variable
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
x = np.linspace(-5,10,1500)
f = stats.uniform.pdf(x,-3,4)
F = stats.uniform.cdf(x,-3,4)
plt.plot(x,f); plt.show()
plt.plot(x,F); plt.show()
# Julia code to generate the PDF and CDF of a uniform random variable
using Distributions
using Plots

u = Uniform(-3, 4)

# pdf
p1 = plot(x -> 0, -5:0.01:1, fillrange=x->pdf(u, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(u, x), -5, 10,
    color=1, linewidth=6, ylims=(0, 0.4), label="PDF")

# cdf
p2 = plot(x -> cdf(u, x), -5:0.01:10,
    linewidth=6, label="CDF")
vline!(p2, [-3, 4],
    linestyle=:dash, color=:green, label=false)

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

CDF of a exponential random variable

% MATLAB code to generate the PDF and CDF of an exponential random variable
pd = makedist('exp',2);

x  = linspace(-5, 10, 1500)';
f  = pdf(pd, x);
x1 = linspace(-5,1,1500)';
f1 = pdf(pd,x1);
figure;
area(x1,f1,'EdgeColor','none','FaceColor',[0.8,0.8,1]);hold on;
plot(x, f, 'LineWidth', 6);
grid on;
axis([-5 10 0 0.6]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);

F = cdf(pd, x);
figure;
plot(x, F, 'LineWidth', 6); hold on;
plot(x(600), F(600), 'ko', 'LineWidth', 2, 'MarkerSize', 10, 'MarkerFaceColor', [0.8,0.8,1]);
grid on;
axis([-5 10 0 1]);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);
# Python code to generate the PDF and CDF
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
x = np.linspace(-5,10,1500)
f = stats.expon.pdf(x,2)
F = stats.expon.cdf(x,2)
plt.plot(x,f); plt.show()
plt.plot(x,F); plt.show()
# Julia code to generate the PDF and CDF of an exponential random variable
using Distributions
using Plots

u = Exponential(2)

# pdf
p1 = plot(x -> 0, -5:0.01:1, fillrange=x->pdf(u, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(u, x), -5, 10,
    color=1, linewidth=6, ylims=(0, 0.6), label="PDF")

# cdf
p2 = plot(x -> cdf(u, x), -5:0.01:10,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

Chapter 4.5

Generate a uniform random variable

% MATLAB code to generate 1000 uniform random numbers
a = 0; b = 1;
X = unifrnd(a,b,[1000,1]);
hist(X);
# Python code to generate 1000 uniform random numbers
import scipy.stats as stats
a = 0; b = 1;
X = stats.uniform.rvs(a,b,size=1000)
plt.hist(X);
# Julia code to generate 1000 uniform random numbers
using Distributions
using Plots

u = Uniform(0, 1) # same as Uniform()
X = rand(u, 1000)
histogram(X)

Mean, variance, median, mode of a uniform random variable

% MATLAB code to compute empirical mean, var, median, mode
X = unifrnd(a,b,[1000,1]);
M = mean(X);
V = var(X);
Med = median(X);
Mod = mode(X);
# Python code to compute empirical mean, var, median, mode
X = stats.uniform.rvs(a,b,size=1000)
M = np.mean(X)
V = np.var(X)
Med = np.median(X)
Mod = stats.mode(X)
# Julia code to compute empirical mean, var, median, mode
using Distributions
u = Uniform(0, 1)
X = rand(u, 1000)

M = mean(X)
V = var(X)
Med = median(X)
Mod = mode(X)

Probability of a uniform random variable

mathbf{P}[0.2 le X le 0.3]
% MATLAB code to compute the probability P(0.2 < X < 0.3)
a = 0; b = 1;
F = unifcdf(0.3,a,b) - unifcdf(0.2,a,b)
# Python code to compute the probability P(0.2 < X < 0.3)
import scipy.stats as stats
a = 0; b = 1;
F = stats.uniform.cdf(0.3,a,b)-stats.uniform.cdf(0.2,a,b)
# Julia code to compute the probability P(0.2 < X < 0.3)
using Distributions

u = Uniform(0, 1)
F = cdf(u, 0.3) - cdf(u, 0.2)

PDF of an exponential random variable

% MATLAB code to generate PDF and CDF of an exponential random variable
lambda1 = 1/2;
lambda2 = 1/5;

x = linspace(0,1,1000);
f1 = pdf('exp',x, lambda1);
f2 = pdf('exp',x, lambda2);

figure(1);
plot(x, f1, '-.', 'LineWidth', 4, 'Color', [0 0.2 0.8]); hold on;
plot(x, f2, 'LineWidth', 4, 'Color', [0.8 0.2 0]);
axis([0 1 0 6]);
legend('\lambda = 2', '\lambda = 5','Location','NE');
grid on;
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);

F1 = cdf('exp',x, lambda1);
F2 = cdf('exp',x, lambda2);
figure(2);
plot(x, F1, '-.', 'LineWidth', 4, 'Color', [0 0.2 0.8]); hold on;
plot(x, F2, 'LineWidth', 4, 'Color', [0.8 0.2 0]);
axis([0 1 0 1.2]);
legend('\lambda = 2', '\lambda = 5','Location','SE');
grid on;
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);
# Python code to generate PDF and CDF of an exponential random variable
lambd1 = 1/2
lambd2 = 1/5
x = np.linspace(0,1,1000)
f1 = stats.expon.pdf(x,scale=lambd1)
f2 = stats.expon.pdf(x,scale=lambd2)
plt.plot(x, f1)
plt.plot(x, f2)

F1 = stats.expon.cdf(x,scale=lambd1)
F2 = stats.expon.cdf(x,scale=lambd2)
plt.plot(x, F1)
plt.plot(x, F2)
# Julia code to generate PDF and CDF of an exponential random variable
using Distributions
using Plots

exp1 = Exponential(1/2)
exp2 = Exponential(1/5)

# Plotting the pdfs
p1 = plot(x -> pdf(exp1, x), 0, 1,
    linewidth=4, linestyle=:dashdot, color=RGB(0,0.2,0.8), label="lambda = 2")
plot!(p1, x -> pdf(exp2, x), 0, 1,
    linewidth=4, color=RGB(0.8,0.2,0), label="lambda = 5")

# Plotting the cdfs
p2 = plot(x -> cdf(exp1, x), 0, 1,
    linewidth=4, linestyle=:dashdot, color=RGB(0,0.2,0.8), label="lambda = 2", legend=:topleft)
plot!(p2, x -> cdf(exp2, x), 0, 1,
    linewidth=4, color=RGB(0.8,0.2,0), label="lambda = 5")

# Figure with both plots
plot(p1, p2, layout=(1, 2),
    size=(600, 400))

Chapter 4.6

PDF and CDF of a Gaussian random variable

% MATLAB code to generate standard Gaussian PDF and CDF
x = linspace(-5,5,1000);
f = normpdf(x,0,1);

x1 = linspace(-5,-1,1000);
f1 = normpdf(x1,0,1);

area(x1,f1,'EdgeColor','none','FaceColor',[0.8,0.8,0.8]);hold on;
plot(x,f,'LineWidth',6,'Color',[0,0,0]);
grid on;
xticks(-5:1:5);
yticks(0:0.1:0.5);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',18);
axis([-5 5 0 0.5]);

figure;
F = normcdf(x,0,1);
plot(x,F,'LineWidth',6,'Color',[0,0,0]); hold on;
plot(x(401), F(401), 'ko', 'LineWidth', 2, 'MarkerSize', 15, 'MarkerFaceColor', [0.8 0.8 0.8]);
grid on;
xticks(-5:1:5);
yticks(0:0.2:1);
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',18);
axis([-5 5 0 1]);
# Python code to generate standard Gaussian PDF and CDF
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
x = np.linspace(-10,10,1000)
f = stats.norm.pdf(x)
F = stats.norm.cdf(x)
plt.plot(x,f); plt.show()
plt.plot(x,F); plt.show()
# Julia code to generate standard Gaussian PDF and CDF
using Distributions
using Plots

normdist = Normal()

p1 = plot(x -> 0, -5:0.01:-1, fillrange=x->pdf(normdist, x),
    linealpha=0, fillcolor=RGB(0.8, 0.8, 1), alpha=0.5, label=false)
plot!(p1, x -> pdf(normdist, x), -5, 5,
    color=1, linewidth=6, ylims=(0, 0.6), label="PDF")

p2 = plot(x -> cdf(normdist, x), -5:0.01:5,
    linewidth=6, label="CDF")

plot(p1, p2, layout=(1, 2),
    size=(600, 300), legend=:topleft)

Skewness and kurtosis of a random variable

% MATLAB code to plot a Gamma distribution
x = linspace(0,30,1000);
theta = 1;
plot(x,gampdf(x,2,theta),'LineWidth',4, 'Color',[0 0 0]); hold on;
plot(x,gampdf(x,5,theta),'LineWidth',4, 'Color',[0.2 0.2 0.2]);
plot(x,gampdf(x,10,theta),'LineWidth',4, 'Color',[0.4 0.4 0.4]);
plot(x,gampdf(x,15,theta),'LineWidth',4, 'Color',[0.6 0.6 0.6]);
plot(x,gampdf(x,20,theta),'LineWidth',4, 'Color',[0.8 0.8 0.8]); hold off;
legend('k = 2', 'k = 5', 'k = 10', 'k = 15', 'k = 20', 'Location', 'NE');
grid on;
set(gcf, 'Position', [100, 100, 600, 300]);
set(gca,'FontWeight','bold','fontsize',14);
% MATLAB code to compute skewness and kurtosis
X = random('gamma',3,5,[10000,1]);
s = skewness(X);
k = kurtosis(X);
# Python code to compute skewness and kurtosis
import scipy.stats as stats
X = stats.gamma.rvs(3,5,size=10000)
s = stats.skew(X)
k = stats.kurtosis(X)
# Julia code to plot a Gamma distribution
using Distributions
using Plots

theta = 1

p1 = plot(size=(600, 300))
plot!(p1, x -> pdf(Gamma(2, theta), x), 0, 30, lw=4, c=RGB(0,0,0), label="k = 2")
plot!(p1, x -> pdf(Gamma(5, theta), x), 0, 30, lw=4, c=RGB(0.2,0.2,0.2), label="k = 5")
plot!(p1, x -> pdf(Gamma(10, theta), x), 0, 30, lw=4, c=RGB(0.4,0.4,0.4), label="k = 10")
plot!(p1, x -> pdf(Gamma(15, theta), x), 0, 30, lw=4, c=RGB(0.6,0.6,0.6), label="k = 15")
plot!(p1, x -> pdf(Gamma(20, theta), x), 0, 30, lw=4, c=RGB(0.8,0.8,0.8), label="k = 20")

# Julia code to compute skewness and kurtosis
X = rand(Gamma(3, 5), 10_000)
s = skewness(X)
k = kurtosis(X)

Chapter 4.8

Generating Gaussians from uniform

Let U be a uniform(0,1) random variable. Apply the transformation

g(U) = F_X^{-1}(U) = sigmaPhi^{-1}(U) + mu.

Then g(U) will be a Gaussian.

% MATLAB code to generate Gaussian from uniform
mu    = 3;
sigma = 2;
U  = rand(10000,1);
gU = sigma*icdf('norm',U,0,1)+mu;
[num,val] = hist(U,30);
figure; bar(val,num,'FaceColor',[0.8, 0.8,0.8]);
set(gcf, 'Position', [100, 100, 600, 400]); grid on;
set(gca,'FontWeight','bold','fontsize',14);

[num,val] = hist(gU,30);
figure; bar(val,num,'FaceColor',[1, 0.8, 0]);
set(gcf, 'Position', [100, 100, 600, 400]); grid on;
set(gca,'FontWeight','bold','fontsize',14);
# Python code to generate Gaussian from uniform
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

mu = 3
sigma = 2
U  = stats.uniform.rvs(0,1,size=10000)
gU = sigma*stats.norm.ppf(U)+mu
plt.hist(U); plt.show()
plt.hist(gU); plt.show()
# Julia code to generate Gaussian from uniform
using Distributions
using Plots

mu = 3
sigma = 2
U = rand(10_000)
gU = sigma * quantile(Normal(), U) .+ mu

p1 = histogram(U, label="U")
p2 = histogram(gU, label="gU")
plot(p1, p2, layout=(1, 2), size=(600, 300), legend=:topleft)