2. Probability Distributions

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from prml.rv import (
    Bernoulli,
    Beta,
    Categorical,
    Dirichlet,
    Gamma,
    Gaussian,
    MultivariateGaussian,
    MultivariateGaussianMixture,
    StudentsT,
    Uniform
)

np.random.seed(1234)

2.1 Binary Variables

In [2]:
model = Bernoulli()
model.fit(np.array([0., 1., 1., 1.]))
print(model)
Bernoulli(
    mu=0.75
)

2.1.1 The beta distributions

In [3]:
x = np.linspace(0, 1, 100)
for i, [a, b] in enumerate([[0.1, 0.1], [1, 1], [2, 3], [8, 4]]):
    plt.subplot(2, 2, i + 1)
    beta = Beta(a, b)
    plt.xlim(0, 1)
    plt.ylim(0, 3)
    plt.plot(x, beta.pdf(x))
    plt.annotate("a={}".format(a), (0.1, 2.5))
    plt.annotate("b={}".format(b), (0.1, 2.1))
plt.show()
In [4]:
beta = Beta(2, 2)
plt.subplot(2, 1, 1)
plt.xlim(0, 1)
plt.ylim(0, 2)
plt.plot(x, beta.pdf(x))
plt.annotate("prior", (0.1, 1.5))

model = Bernoulli(mu=beta)
model.fit(np.array([1]))
plt.subplot(2, 1, 2)
plt.xlim(0, 1)
plt.ylim(0, 2)
plt.plot(x, model.mu.pdf(x))
plt.annotate("posterior", (0.1, 1.5))

plt.show()
In [5]:
print("Maximum likelihood estimation")
model = Bernoulli()
model.fit(np.array([1]))
print("{} out of 10000 is 1".format(model.draw(10000).sum()))

print("Bayesian estimation")
model = Bernoulli(mu=Beta(1, 1))
model.fit(np.array([1]))
print("{} out of 10000 is 1".format(model.draw(10000).sum()))
Maximum likelihood estimation
10000 out of 10000 is 1
Bayesian estimation
6649 out of 10000 is 1

2.2 Multinomial Variables

In [6]:
model = Categorical()
model.fit(np.array([[0, 1, 0], [1, 0, 0], [1, 0, 0], [0, 0, 1]]))
print(model)
Categorical(
    mu=[0.5  0.25 0.25]
)

2.2.1 The Dirichlet distribution

In [7]:
mu = Dirichlet(alpha=np.ones(3))
model = Categorical(mu=mu)
print(model)

model.fit(np.array([[1., 0., 0.], [1., 0., 0.], [0., 1., 0.]]))
print(model)
Categorical(
    mu=Dirichlet(
        alpha=[1. 1. 1.]
    )
)
Categorical(
    mu=Dirichlet(
        alpha=[3. 2. 1.]
    )
)

2.3 The Gaussian Distribution

In [8]:
uniform = Uniform(low=0, high=1)
plt.figure(figsize=(10, 5))

plt.subplot(1, 3, 1)
plt.xlim(0, 1)
plt.ylim(0, 5)
plt.annotate("N=1", (0.1, 4.5))
plt.hist(uniform.draw(100000), bins=20, density=True)

plt.subplot(1, 3, 2)
plt.xlim(0, 1)
plt.ylim(0, 5)
plt.annotate("N=2", (0.1, 4.5))
plt.hist(0.5 * (uniform.draw(100000) + uniform.draw(100000)), bins=20, density=True)

plt.subplot(1, 3, 3)
plt.xlim(0, 1)
plt.ylim(0, 5)
sample = 0
for _ in range(10):
    sample = sample + uniform.draw(100000)
plt.annotate("N=10", (0.1, 4.5))
plt.hist(sample * 0.1, bins=20, density=True)

plt.show()

2.3.4 Maximum Likelihood for the Gaussian

In [9]:
X = np.random.normal(loc=1., scale=2., size=(100, 2))
gaussian = MultivariateGaussian()
gaussian.fit(X)
print(gaussian)

x, y = np.meshgrid(
    np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
p = gaussian.pdf(
    np.array([x, y]).reshape(2, -1).T).reshape(100, 100)
plt.scatter(X[:, 0], X[:, 1], facecolor="none", edgecolor="steelblue")
plt.contour(x, y, p)
plt.show()
MultivariateGaussian(
    mu=[0.91852581 1.17919155]
    cov=[[4.29224408 0.1551223 ]
 [0.1551223  3.58170912]]
)

2.3.6 Bayesian inference for the Gaussian

In [10]:
mu = Gaussian(0, 0.1)
model = Gaussian(mu, 0.1)

x = np.linspace(-1, 1, 100)
plt.plot(x, model.mu.pdf(x), label="N=0")

model.fit(np.random.normal(loc=0.8, scale=0.1, size=1))
plt.plot(x, model.mu.pdf(x), label="N=1")

model.fit(np.random.normal(loc=0.8, scale=0.1, size=1))
plt.plot(x, model.mu.pdf(x), label="N=2")

model.fit(np.random.normal(loc=0.8, scale=0.1, size=8))
plt.plot(x, model.mu.pdf(x), label="N=10")

plt.xlim(-1, 1)
plt.ylim(0, 5)
plt.legend()
plt.show()
In [11]:
x = np.linspace(0, 2, 100)
for i, [a, b] in enumerate([[0.1, 0.1], [1, 1], [2, 3], [4, 6]]):
    plt.subplot(2, 2, i + 1)
    gamma = Gamma(a, b)
    plt.xlim(0, 2)
    plt.ylim(0, 2)
    plt.plot(x, gamma.pdf(x))
    plt.annotate("a={}".format(a), (1, 1.6))
    plt.annotate("b={}".format(b), (1, 1.3))
plt.show()
In [12]:
tau = Gamma(a=1, b=1)
model = Gaussian(mu=0, tau=tau)
print(model)

model.fit(np.random.normal(scale=1.414, size=100))
print(model)
Gaussian(
    mu=0
    tau=Gamma(
        a=1
        b=1
    )
    var=None
)
Gaussian(
    mu=0
    tau=Gamma(
        a=51.0
        b=94.73272357871433
    )
    var=None
)

2.3.7 Student's t-distribution

In [13]:
X = np.random.normal(size=20)
X = np.concatenate([X, np.random.normal(loc=20., size=3)])
plt.hist(X.ravel(), bins=50, density=1., label="samples")

students_t = StudentsT()
gaussian = Gaussian()

gaussian.fit(X)
students_t.fit(X)

print(gaussian)
print(students_t)

x = np.linspace(-5, 25, 1000)
plt.plot(x, students_t.pdf(x), label="student's t", linewidth=2)
plt.plot(x, gaussian.pdf(x), label="gaussian", linewidth=2)
plt.legend()
plt.show()
Gaussian(
    mu=2.2695020512383577
    var=46.95748684941486
    tau=0.02129585859666723
)
StudentsT(
    mu=-0.1946342109447806
    tau=1.4665461746068318
    dof=0.9028354354811657
)

2.3.9 Mixture of Gaussians

In [14]:
x1 = np.random.normal(size=(100, 2))
x1 += np.array([-5, -5])
x2 = np.random.normal(size=(100, 2))
x2 += np.array([5, -5])
x3 = np.random.normal(size=(100, 2))
x3 += np.array([0, 5])
X = np.vstack((x1, x2, x3))

model = MultivariateGaussianMixture(n_components=3)
model.fit(X)
print(model)

x_test, y_test = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
X_test = np.array([x_test, y_test]).reshape(2, -1).transpose()
probs = model.pdf(X_test)
Probs = probs.reshape(100, 100)
plt.scatter(X[:, 0], X[:, 1])
plt.contour(x_test, y_test, Probs)
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()
MultivariateGaussianMixture(
    mu=[[ 5.06087392 -5.07813706]
 [-4.9724865  -5.09928156]
 [-0.05584106  4.99523381]]
    cov=[[[ 1.11567912 -0.01717074]
  [-0.01717074  1.04457722]]

 [[ 0.89251355 -0.0138146 ]
  [-0.0138146   1.07986159]]

 [[ 0.81797404  0.03778106]
  [ 0.03778106  0.93690783]]]
    coef=[0.33333333 0.33333333 0.33333333]
)
In [ ]: