Python, MATLAB, Julia, R code: Chapter 3

% 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)
plt.show()

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

% 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)
plt.show()

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

% 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);
plt.show()

# 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
X <- rexp(k, rate=lambda)
hist(X, 200)

% 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 <- rexp(n, rate=lambda)
m <- 6:200
J <- numeric(195)
for (i in 1:195) {
    num <- hist(X, breaks=seq(0, max(X), l=m[i]), plot=FALSE)$counts
    h <- n/m[i]
    J[i] = 2/((n-1)*h)-((n+1)/((n-1)*h))*sum((num/n)^2)
}
plot(m, J, type="l", lwd=3)

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

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

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

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

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

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

% 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');

# Python code to generate Erdos Renyi Graph
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
A = np.random.rand(40, 40) < 0.3
A = np.triu(A, 1)
A = A + A.T
G = nx.convert_matrix.from_numpy_matrix(A)
fig, ax = plt.subplots()
nx.draw(G, with_labels=True, ax=ax)
limits = plt.axis("on")
ax.tick_params(left=True, bottom=True, labelleft=True, labelbottom=True)
plt.show()

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

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

% 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")

% 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);

# Python code to plot the Poisson PMF
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
lambda_set = [1, 4, 10]
p = np.zeros((20, 3))
k = np.arange(0, 20)
for i in range(0, 3):
    lambd = lambda_set[i]
    p[:,i] = lambd**k/factorial(k)*np.exp(-lambd)

plt.plot(k, p[:,0], 'bo')
plt.vlines(k, 0 , p[:,0], colors='b', lw=2)
plt.plot(k, p[:,1], 'go')
plt.vlines(k, 0 , p[:,1], colors='g', lw=2)
plt.plot(k, p[:,2], 'yo')
plt.vlines(k, 0 , p[:,2], colors='y', lw=2)
labels = ["lambda = 1", "lambda = 4", "lambda = 10"]
plt.legend(labels=labels)
plt.grid()
plt.show()

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

# R code to plot the Poisson PMF
library(pracma)
lambda_set = c(1, 4, 10)
p = zeros(21, 3)
k = seq(0, 20)
for (i in 1:3) {
    p[,i] = dpois(k, lambda_set[i])
}

plot(k, p[,1], type='h', lwd=3, col="blue")
grid()
legend(15, 0.36, legend=c("lambda  = 1", "lambda  = 4", "lambda  = 10"), col=c("blue", "darkgreen", "orange"), lwd=3)
lines(k, p[,1], type='p', cex=2, pch=16, col="blue")
lines(k, p[,2], type='h', lwd=3, col="darkgreen")
lines(k, p[,2], type='p', cex=2, pch=16, col="darkgreen")
lines(k, p[,3], type='h', lwd=3, col="orange")
lines(k, p[,3], type='p', cex=2, pch=16, col="orange")

% 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);

# Python code to plot the Poisson CDF
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
lambda_set = [1, 4, 10]
p = np.zeros((20, 3))
k = np.arange(0, 20)
for i in range(0, 3):
    lambd = lambda_set[i]
    p[:,i] = lambd**k/factorial(k)*np.exp(-lambd)

plt.step(k, np.cumsum(p[:,0]), 'b')
plt.step(k, np.cumsum(p[:,1]), 'g')
plt.step(k, np.cumsum(p[:,2]), 'y')
plt.plot(k, np.cumsum(p[:,0]), 'bo')
plt.plot(k, np.cumsum(p[:,1]), 'go')
plt.plot(k, np.cumsum(p[:,2]), 'yo')
labels = ["lambda = 1", "lambda = 4", "lambda = 10"]
plt.legend(labels=labels)
plt.grid()
plt.show()

# 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(lambda _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)

# R code to plot the Poisson CDF
library(pracma)
lambda_set = c(1, 4, 10)
p = zeros(21, 3)
k = seq(0, 20)
for (i in 1:3) {
    p[,i] = ppois(k, lambda_set[i])
}

plot(k, p[,3], type='s', lwd=3, col="orange")
grid()
legend(15, 0.36, legend=c("lambda  = 1", "lambda  = 4", "lambda  = 10"), col=c("blue", "darkgreen", "orange"), lwd=3)
lines(k, p[,3], type='p', cex=2, pch=16, col="orange")
lines(k, p[,2], type='s', lwd=3, col="darkgreen")
lines(k, p[,2], type='p', cex=2, pch=16, col="darkgreen")
lines(k, p[,1], type='s', lwd=3, col="blue")
lines(k, p[,1], type='p', cex=2, pch=16, col="blue")

% 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
lambda = n*p
x = 0:120
y = dbinom(x, n, p)
z = dpois(x, lambda)

plot(x, y, lwd=2, type="h", col="darkolivegreen3", xlab="k", ylab="Probability")
legend(63, 0.053, legend=c("Binomial, n=5000, p=0.01", "Poisson, lambda = 50"), col=c("darkolivegreen3", "blue"), lwd=3)
grid()
lines(x, y, lwd=2, cex=1, type="p", col="darkolivegreen3")
lines(x, z, lwd=2, type="l", col="blue")

% 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
x = cv2.imread('cameraman.tif', 0)/255
width = len(x)

alpha = 10
lam = alpha * x
X1 = np.random.poisson(lam=lam, size=(width, width))

alpha = 100
lam = alpha * x
X2 = np.random.poisson(lam=lam, size=(width, width))

alpha = 1000
lam = alpha * x
X3 = np.random.poisson(lam=lam, size=(width, width))

plt.figure(1)
plt.imshow(X1, cmap='gray')
plt.title("alpha = 10")
plt.figure(2)
plt.imshow(X2, cmap='gray')
plt.title("alpha = 100")
plt.figure(3)
plt.imshow(X3, cmap='gray')
plt.title("alpha = 1000")
plt.show()

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

# R code to demonstrate the photon shot noise
library(imager)
x = as.data.frame(load.image("cameraman.tif"))
x = xtabs(value ~ x+y, data=x)
width = sqrt(length(x))

alpha = 10
lambda = alpha * x
X1 = rpois(length(lambda), lambda)
X1 = array(X1, c(width, width))

alpha = 100
lambda = alpha * x
X2 = rpois(length(lambda), lambda)
X2 = array(X2, c(width, width))

alpha = 1000
lambda = alpha * x
X3 = rpois(length(lambda), lambda)
X3 = array(X3, c(width, width))

# Flip matrix since `image()` reads the matrix bottom up
flip_matrix = function(m) m[,nrow(m):1]

# Plot images together
layout(matrix(1:3, 1, 3), respect=TRUE)
image(flip_matrix(X1), col=gray.colors(255), main = expression(paste(alpha, " = 10")))
image(flip_matrix(X2), col=gray.colors(255), main = expression(paste(alpha, " = 100")))
image(flip_matrix(X3), col=gray.colors(255), main = expression(paste(alpha, " = 1000")))