Python, MATLAB, Julia, R code: Chapter 4

# 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)

# R code to generate the PDF and CDF
library(EnvStats)

x = seq(-5, 10, (5+10)/1000)
x = sort(x)
f = dnormMix(x, 0, 1, 5, 1, p.mix=c(0.7))
plot(x,f, xlim=c(-5,10), ylim=c(0,0.4), type="n")
lines(x, f, lwd=4, col="blue")
polygon(c(min(x), x[x<=0.8]), c(f[x<=0.8], 0), col="lightblue")
grid()

F = pnormMix(x, 0, 1, 5, 1, p.mix=c(0.7))
plot(x, F, type="n")
lines(x, F, lwd=4, col="blue")
grid()

# 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)

# R code to generate the PDF and CDF
x = seq(-5, 10, (5+10)/1500)
f = dunif(x, -3, 4)
F = punif(x, -3, 4)
plot(x, f, type="n")
lines(x, f, lwd=5)
grid()
plot(x, F, type="n")
lines(x, F, lwd=5)
grid()

# 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)

# R code to generate the PDF and CDF
x = seq(-5, 10, (5+10)/1500)
f = dexp(x, 1/2)
F = pexp(x, 1/2)
# PDF
plot(x, f, type="n")
lines(x, f, lwd=5, col="blue")
grid()

# CDF
plot(x, F, type="n")
lines(x, F, lwd=5, col="blue")
grid()

# 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)

# R code to generate 1000 uniform random numbers
a = 0; b = 1;
X = runif(1000, a, b)
hist(X)
grid()

# 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)

# R code to computer empirical mean, var, median, mode
library(pracma)
a = 0; b = 1;
X = runif(1000, a, b)
M = mean(X)
V = var(X)
Med = median(X)
Mod = Mode(X)

# 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)

# R code to compute the probability P(0.2 < X < 0.3)
a = 0; b = 1;
F = punif(0.3, a, b) - punif(0.2, a, b)

# 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))

# R code to generate the PDF and CDF of an exponential random variable
lambda1 = 1/2
lambda2 = 1/5
x = seq(0, 1, (0+1)/1000)
f1 = dexp(x, 1/lambda1)
f2 = dexp(x, 1/lambda2)
plot(x, f2, type="n")
lines(x, f1, lwd=4, col="blue")
lines(x, f2, lwd=4, col="red")
legend(0, 1, legend=c(expression(paste(lambda, "=5")), expression(paste(lambda, "=2"))), col=c("red", "blue"), lty=1:1)
grid()

F1 = pexp(x, 1/lambda1)
F2 = pexp(x, 1/lambda2)
plot(x, F2, type="n")
lines(x, F1, lwd=4, col="blue")
lines(x, F2, lwd=4, col="red")
legend(0, 1, legend=c(expression(paste(lambda, "=5")), expression(paste(lambda, "=2"))), col=c("red", "blue"), lty=1:1)
grid()

# 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)

# R code to generate a Gaussian PDF
x = seq(-10, 10, (10+10)/1000)
mu = 0; sigma = 1;
f = dnorm(x, mu, sigma)
plot(x, f, type="n")
lines(x, f, lwd=4, col="blue")
grid()

# R code to generate standard Gaussian PDF and CDF
x = seq(-5, 5, (5+5)/1000)
f = dnorm(x)
F = pnorm(x)
plot(x, f, type = "n")
lines(x, f, lwd=6)
grid()
plot(x, F, type = "n")
lines(x, F, lwd=6)
grid()

# 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)

# R code to plot a Gamma distribution
x = seq(0, 30, (0+30)/1000)
theta = 1
plot(x, dgamma(x, 2, theta), type = "n")
lines(x, dgamma(x, 2, theta), lwd = 4)
lines(x, dgamma(x, 5, theta), lwd = 4, col = "#333333")
lines(x, dgamma(x, 10, theta), lwd = 4, col = "#666666")
lines(x, dgamma(x, 15, theta), lwd = 4, col = "#999999")
lines(x, dgamma(x, 20, theta), lwd = 4, col = "#CCCCCC")
legend(23, 0.36, legend=c("k = 2", "k = 5", "k = 10", "k = 15", "k = 20"), col=c("1", "#333333", "#666666", "#999999", "#CCCCCC"), lwd = 4,lty=1:1)
grid()

library(e1071)
X = rgamma(10000, 3, 5)
s = skewness(X)
k = kurtosis(X)

# 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)

# R code to generate Gaussian from uniform
library(HDInterval)
mu = 3
sigma = 2
U = runif(10000, 0, 1)
gU = sigma * inverseCDF(U, pnorm) + mu;
hist(U)
grid()

hist(gU, col="yellow")
grid()