Python, MATLAB, Julia, R code: Chapter 3

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

Acknowledgement: The R code is written by contributors listed here.

Chapter 3.2

Histogram of the English alphabets

Data download: ch3_data_english.txt (1KB)

% MATLAB code to generate the histogram
load('ch3_data_English');
bar(f/100,'FaceColor',[0.9,0.6,0.0]);
xticklabels({'a','b','c','d','e','f','g','h','i','j','k','l',...
    'm','n','o','p','q','r','s','t','u','v','w','x','y','z'});
xticks(1:26);
yticks(0:0.02:0.2);
axis([1 26 0 0.13]);
# Python code generate the histogram
import numpy as np
import matplotlib.pyplot as plt
f = np.loadtxt('./ch3_data_english.txt')
n = np.arange(26)
plt.bar(n, f/100)
ntag = ['a','b','c','d','e','f','g','h','i','j','k','l','m',...
    'n','o','p','q','r','s','t','u','v','w','x','y','z']
plt.xticks(n, ntag)
# Julia code generate the histogram
using Plots,DelimitedFiles
f = readdlm("ch3_data_english.txt")
bar(f/100,xticks = (1:26,'a':'z'),yticks = 0.0:0.02:0.1,label = false)
# R code generate the histogram

f <- read.table('./ch3_data_english.txt')
f <- f/100
n <- c(1:26)
ntag <- c('a','b','c','d','e','f','g','h','i','j','k','l','m',
    'n','o','p','q','r','s','t','u','v','w','x','y','z')

f$index <- n
f$labels <- ntag

barplot(f$V1~f$index, xaxt = 'n', xlab = "Letters", ylab = "Values")
axis(1, at=1:26, labels=ntag)

Histogram of the throwing a die

% MATLAB code to generate the histogram
x = [1 2 3 4 5 6];
q = randi(6,100,1);
[num,val] = hist(q,x-0.5);
bar(num/100,'FaceColor',[0.8, 0.8,0.8]);
axis([0 7 0 0.24]);
# Python code generate the histogram
import numpy as np
import matplotlib.pyplot as plt
q = np.random.randint(7,size=100)
plt.hist(q+0.5,bins=6)
# Julia code generate the histogram
using Plots,Random
q = rand(1:6,100)
histogram(q,bins = 6,normalize = true,label = "N = 100")
# R code generate the histogram
q <- sample(1:6, 100, replace = T)
hist(q + 0.5, 6)

Histogram of an exponential random variable

% MATLAB code used to generate the plots
lambda = 1;
k      = 1000;
X      = exprnd(1/lambda,[k,1]);
[num,val] = hist(X,200);
bar(val,num,'FaceColor',[1, 0.5,0.5]);
# Python code used to generate the plots
import numpy as np
import matplotlib.pyplot as plt
lambd = 1
k     = 1000
X     = np.random.exponential(1/lambd, size=k)
plt.hist(X,bins=200);
# Julia code used to generate the plots
using Plots,Random,Distributions
lambda = 1
X = rand(Exponential(1/lambda),1_000)
histogram(X,bins = 200,label = "K = 200")
# R code used to generate the plots
lambda <- 1
k <- 1000
set = runif(n=k, min=0, max=1)
freq = -1/lambda * log(1-set)
hist(freq, 200)

Cross validation loss

% MATLAB code to perform the cross validation
lambda = 1;
n = 1000;
X = exprnd(1/lambda,[n,1]);
m = 6:200;
J = zeros(1,195);

for i=1:195
    [num,binc] = hist(X,m(i));
    h = n/m(i);
    J(i) = 2/((n-1)*h)-((n+1)/((n-1)*h))*sum( (num/n).^2 );
end

plot(m,J,'LineWidth',4,'Color',[0.9,0.2,0.0]);
# Python code to perform the cross validation
import numpy as np
import matplotlib.pyplot as plt
lambd = 1
n     = 1000
X     = np.random.exponential(1/lambd, size=n)
m     = np.arange(5,200)
J     = np.zeros((195))

for i in range(0,195):
  hist,bins = np.histogram(X,bins=m[i])
  h = n/m[i]
  J[i] = 2/((n-1)*h)-((n+1)/((n-1)*h))*np.sum((hist/n)**2)

plt.plot(m,J);
# Julia code to perform the cross validation

using Plots,Random,Distributions,StatsBase
lambda = 1
n = 1_000
X = rand(Exponential(1/lambda),n)
m = 6:200
J = zeros(195)
for i = 1:195
    hist = fit(Histogram,X,range(minimum(X),maximum(X),length = m[i]+1))
    h = n/m[i]
    J[i] = 2/((n-1)*h)-((n+1)/((n-1)*h))*sum( (hist.weights/n).^2 )
end

plot(m,J,label = false)
# R code to perform the cross validation
lambda <- 1
n <- 1000
X <- runif(n=n, min=0, max=1)
freq = -1/lambda * log(1-X)
m <- c(5:200)
J <- replicate(195, 0)

for (i in 1:195) {
    h <- n/m[i]
    J[i] = (2/((n-1)*h)-((n+1)/((n-1)*h)))*sum((hist(freq, m[i])/n)^2)
}
plot(J, m)

Chapter 3.4

Mean of a vector

Expectation mathbf{E}[X] where X is uniform from [0,1]

% MATLAB code to compute the mean of a dataset
X  = rand(10000,1);
mX = mean(X);
# Python code to compute the mean of a dataset
import numpy as np
X  = np.random.rand(10000)
mX = np.mean(X)
# Julia code to compute the mean of a dataset
using Random,Statistics
X = rand(10_000)
mean(X)
# R code to compute the mean of a dataset
X <- runif(n=10000, min=0, max=1)
print(mean(X))

Mean from PMF

% MATLAB code to compute the expectation
p = [0.25 0.5 0.25];
x = [0 1 2];
EX = sum(p.*x);
# Python code to compute the expectation
import numpy as np
p = np.array([0.25, 0.5, 0.25])
x = np.array([0, 1, 2])
EX = np.sum(p*x)
# Julia code to compute the expectation
p = [0.25,0.5,0.25]
x = [0,1,2]
EX = sum(p .* x)
# R code to compute the expectation
p = c(0.25, 0.5, 0.25)
x = c(0, 1, 2)
EX = sum(p*x)
print(EX)

Mean of geometric random variable

The expectation of a random variable X where

p_X(k) = frac{1}{2^k} quad k = 1,2,3ldots.
% MATLAB code to compute the expectation
k = 1:100;
p = 0.5.^k;
EX = sum(p.*k);
# Python code to compute the expectation
import numpy as np
k = np.arange(100)
p = np.power(0.5,k)
EX = np.sum(p*k)
# Julia code to compute the expectation
k = 1:100
p = 0.5.^k
EX = sum(p .* k)
# R code to compute the expectation
k = c(1:100)
p = 0.5 ^ k
EX = sum(p*k)
print(EX)

Chapter 3.5 Common Discrete Random Variables

Bernoulli Random Variables

% MATLAB code to generate 1000 Bernoulli random variables
p = 0.5;
n = 1;
X = binornd(n,p,[1000,1]);
[num, ~] = hist(X, 10);
bar(linspace(0,1,10), num,'FaceColor',[0.4, 0.4, 0.8]);
# Python code to generate 1000 Bernoulli random variables
import numpy as np
import matplotlib.pyplot as plt
p = 0.5
n = 1
X = np.random.binomial(n,p,size=1000)
plt.hist(X,bins='auto')
# Julia code to generate 1000 Bernoulli random variables
using Plots,Distributions
p = 0.5
n = 1
X = rand(Binomial(n,p),1_000)
histogram(X,bins = 2,label = false)
# R code to generate 1000 Bernoulli random variables
p <- 0.5
n <- 1
X <- rbinom(1000, n, p)
hist(X)

Alternatively, we can call the scipy.stats library

# Python code to call scipy.stats library
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
p = 0.5
X  = stats.bernoulli.rvs(p,size=1000)
plt.hist(X,bins='auto');

To generate Bernoulli random variable objects, we can call the scipy.stats library

# Python code to generate a Bernoulli rv object
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
p = 0.5
rv = stats.bernoulli(p)
mean, var = rv.stats(moments='mv')
print(mean, var)

Erdos Renyi Graph

% MATLAB code to generate Erdos Renyi Graph
A = rand(40,40)<0.3;
A = triu(A,1);
A = A+A';
figure(1);
G = graph(A);
plot(G, 'LineWidth', 2, 'EdgeColor', [0.1 0.6 0], 'MarkerSize',10);
set(gcf, 'Position', [100, 100, 300, 300]);
title('p = 0.3');
# Julia code to generate Erdos Renyi Graph
using Plots,Random,LinearAlgebra,GraphRecipes
A = rand(40,40) .< 0.1
A = triu(A,1)
A = A + A'
graphplot(A)

Binomial Random Variables

% MATLAB code to generate 5000 Binomial random variables
p = 0.5;
n = 10;
X = binornd(n,p,[5000,1]);
[num, ~] = hist(X, 10);
bar( num,'FaceColor',[0.4, 0.4, 0.8]);
# Python code to generate 5000 Binomial random variables
import numpy as np
import matplotlib.pyplot as plt
p = 0.5
n = 10
X = np.random.binomial(n,p,size=5000)
plt.hist(X,bins='auto');
# Julia code to generate 5000 Binomial random variables
using Plots,Distributions
p = 0.5
n = 10
X = rand(Binomial(n,p),5_000)
histogram(X,label = false)
# R code to generate 5000 Binomial random variables
p <- 0.5
n <- 10
X <- rbinom(5000, n, p)
hist(X)

Binomial CDF

% MATLAB code to plot CDF of a binomial random variable
x = 0:10;
p = 0.5;
n = 10;
F = binocdf(x,n,p);
figure; stairs(x,F,'.-','LineWidth',4,'MarkerSize',30);
# Python code to plot CDF of a binomial random variable
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
p = 0.5
n = 10
rv = stats.binom(n,p)
x  = np.arange(11)
F  = rv.cdf(x)
plt.plot(x, F, 'bo', ms=10);
plt.vlines(x, 0, F, colors='b', lw=5, alpha=0.5);
# Julia code to plot CDF of a binomial random variable
using Plots,Distributions
x = 0:10
p = 0.5
n = 10
F = cdf.(Binomial(n,p),x)
plot(x,F,linetype = :steppost,label = false,markershape = :circle)
# R code to plot CDF of a binomial random variable
p <- 0.5
n <- 10
x <- 0:n
plot(x, pbinom(x, size = n, prob = p), type="h")

Poisson CDF

% MATLAB code to plot the Poisson PMF
lambda_set = [1 4 10];
for i=1:3
    lambda = lambda_set(i);
    k = 0:1:20;
    p(:,i) = lambda.^k./factorial(k).*exp(-lambda);
end
figure;
stem(k, p(:,1),'b.','LineWidth',2,'MarkerSize',30); hold on;
stem(k, p(:,2),'.','LineWidth',2,'MarkerSize',30, 'Color', [0.1 0.5 0]);
stem(k, p(:,3),'.','LineWidth',2,'MarkerSize',30, 'Color', [1 0.5 0]);
legend('\lambda = 1', '\lambda = 4', '\lambda = 10');
grid on;
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);
# Julia Code for the Poisson PMF
using Plots,Distributions
lambda_set = [1,4,10]
p = zeros(21,3)
k = 0:1:20
for i = 1:3
    p[:,i] = pdf.(Poisson(lambda_set[i]),k)
end
plot(k,p[:,1],line = :stem,markershape = :circle,label = "lambda = 1")
plot!(k,p[:,2],line = :stem,markershape = :circle,label = "lambda = 4")
plot!(k,p[:,3],line = :stem,markershape = :circle,label = "lambda = 10")
% MATLAB code to plot the Poisson CDF
lambda_set = [1 4 10];
for i=1:3
    lambda = lambda_set(i);
    k = 0:1:20;
    p(:,i) = lambda.^k./factorial(k).*exp(-lambda);
end
figure;
stairs(k, cumsum(p(:,1)), 'b.-','MarkerSize',30, 'LineWidth', 2); hold on;
stairs(k, cumsum(p(:,2)), '.-','MarkerSize',30, 'Color', [0.1 0.5 0], 'LineWidth', 2);
stairs(k, cumsum(p(:,3)), '.-','MarkerSize',30, 'Color', [1 0.5 0], 'LineWidth', 2);
legend('\lambda = 1', '\lambda = 4', '\lambda = 10', 'Location', 'SE');
grid on;
set(gcf, 'Position', [100, 100, 600, 400]);
set(gca,'FontWeight','bold','fontsize',14);
# Julia Code for Poisson CDF
using Plots,Distributions
lambda_set = [1,4,10]
p = zeros(21,3)
k = 0:1:20
for i = 1:3
    p[:,i] = cdf.(Poisson(λ_set[i]),k)
end
plot(k,p[:,1],line = :steppost,markershape = :circle,label = "lambda = 1")
plot!(k,p[:,2],line = :steppost,markershape = :circle,label = "lambda = 4")
plot!(k,p[:,3],line = :steppost,markershape = :circle,label = "lambda = 10",legend = :bottomright)

Poisson-Binomial approximation

% MATLAB code to approximate binomial using Poisson
n = 5000;
p = 0.01;
lambda = n*p;
x = 0:120;
y = binopdf(x,n,p);
z = poisspdf(x,lambda);
stem(x,y,'Color',[0.5,0.7,0],'LineWidth', 2); hold on;
plot(x,z,'b','LineWidth',2);
legend('Binomial, n = 5000, p = 0.01','Poisson, \lambda = 50');
xlabel('k');
ylabel('Probability');
grid on;
set(gcf, 'Position', [100, 100, 800, 300]);
set(gca,'FontWeight','bold','fontsize',14);
# Python code to approximate binomial using Poisson
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
n = 5000; p = 0.01
rv1 = stats.binom(n,p)
X   = rv1.rvs(size=10000)
plt.figure(1); plt.hist(X,bins=np.arange(0,100));
rv2 = stats.poisson(n*p)
f   = rv2.pmf(bin)
plt.figure(2); plt.plot(f);
# Julia code to approximate binomial using Poisson
using Plots,Distributions
n = 5_000
p = 0.01
lambda = n*p
x = 0:120
y = pdf.(Binomial(n,p),x)
z = pdf.(Poisson(lambda),x)
plot(x,y,line = :stem,markershape = :circle,label = "Binomial,n = 5000,p = 0.01")
plot!(x,z,label = "Poisson, lambda = 50")
# R code to approximate binomial using Poisson
n <- 5000
p <- 0.01
x <- rbinom(10000, 5000, 0.01)
pois <- ppois(n*p, x)
plot(x, pois)

Photon Shot Noise

Image download: cameraman.tif (65KB)

% MATLAB to demonstrate the photon shot noise
x = im2double(imread('cameraman.tif'));

alpha = 10;
lambda = alpha*x;
X1 = poissrnd(lambda);

alpha = 100;
lambda = alpha*x;
X2 = poissrnd(lambda);

alpha = 1000;
lambda = alpha*x;
X3 = poissrnd(lambda);

figure;
subplot(131); imshow(X1,[]); title('\alpha = 10');
subplot(132); imshow(X2,[]); title('\alpha = 100');
subplot(133); imshow(X3,[]); title('\alpha = 1000');
# Python to demonstrate the photon shot noise
import numpy as np
import matplotlib.pyplot as plt
import cv2
x0 = cv2.imread('./cameraman.tif', 0)/255
plt.figure(1); plt.imshow(x0,cmap='gray');
X  = np.random.poisson(10*x0)
plt.figure(2); plt.imshow(X, cmap='gray');
# Julia code to demonstrate the photon shot noise
using Distributions,Images,ImageView,Random
x = Float64.(load("cameraman.tif"))

alpha = 10
lambda = alpha*x
X1 = rand.(Poisson.(lambda))
imshow(X1)

alpha = 100
lambda = alpha*x
X1 = rand.(Poisson.(lambda))
imshow(X1)

alpha = 1_000
lambda = alpha*x
X1 = rand.(Poisson.(lambda))
imshow(X1)