The alpha that I create is not created in a crazily efficient way I know, but the main bottleneck in this code is in what I call the "Maximization step" (see comments in code). There are double for-loops twice (with list comprehensions inside), and the reason for this is that I need to use the computed vector from the first double for-loop in the second double for-loop. Here's some proof of my claim, I tested how long the different parts of the code took;
Finishing the vector alphas took 0.1451436 s.
Finishing the first bottleneck for loops took 13.7402911 s
Finishing the second bottleneck for loops took 13.1357222 s
A single complete loop in the while step took 26.9219609 s
There should be a more numpy-focused way to update the three vectors mu, pi, sigma_squared, I just can't seem to figure out how. Any help to make this quicker would be greatly appreciated.
import numpy as np
import random as rnd
from scipy.stats import norm
no_distributions = 2 # range of L
range_of_i = 255
range_of_m = 100
iterations = 100
mu = [5.7e-04, 5.9e-04]
sigma_squared = [9.5e-07, 9.5e-07]
weights = [0.49, 0.51]
alphas = []
for i in range(range_of_i):
alpha = []
for m in range(range_of_m):
sub_populations = [rnd.gauss(mu[i], np.sqrt(sigma_squared[i])) for i in range(no_distributions)]
alpha.append(rnd.choices(sub_populations, weights=weights)[0]) # index 0 to add float instead of list
alphas.append(alpha)
threshold = 10
thresh_limit = 1e-4
initial_parameters = np.array([[-1, 1], [0.2, 0.2], [0.5, 0.5]])
# we iterate between the expectation and maximization steps until convergence
while threshold > thresh_limit:
# mu, pi and sigma_squared need to be numpy arrays because we want to add vectors together the numpy way
mu_initial = initial_parameters[0, :]
sigma_squared_initial = initial_parameters[1, :]
pi_initial = initial_parameters[2, :]
mu_updated = np.zeros(no_distributions)
pi_updated = np.zeros(no_distributions)
sigma_squared_updated = np.zeros(no_distributions)
indicator_normalized = []
# Expectation step
for i in range(range_of_i):
indicator = [pi_initial[l] * norm.pdf(alphas[i], mu_initial[l], np.sqrt(sigma_squared_initial[l])) for l
in range(no_distributions)]
indicator_normalized.append([l / sum(indicator) for l in indicator])
indicator_normalized = np.array(indicator_normalized)
indicator_normalized = indicator_normalized.transpose(0, 2, 1)
# Maximization step
for i in range(range_of_i):
for j in range(np.array(alphas).shape[1]):
mu_updated += np.array(
[indicator_normalized[i][j][l] * alphas[i][j] / sum(sum(indicator_normalized))[l] for l in
range(no_distributions)])
pi_updated += np.array(
[indicator_normalized[i][j][l] / (range_of_i * range_of_m)
for l in range(no_distributions)])
# same for loop again needed because we want to use the complete mu_vector to calculate sigma_squared
for i in range(range_of_i):
for j in range(np.array(alphas).shape[1]):
sigma_squared_updated += np.array([indicator_normalized[i][j][l] * (
alphas[i][j] - mu_updated[l]) ** 2 / sum(sum(indicator_normalized))[l] for l in
range(no_distributions)])
# update parameters and calculate the new threshold to see if convergence has been met
parameters_updated = np.array([mu_updated, sigma_squared_updated, pi_updated])
threshold = sum(sum(np.abs(parameters_updated - initial_parameters)))
initial_parameters = parameters_updated
print (parameters_updated)
The code basically implements an EM-algorithm where the purpose is to find better estimates (better yet, the MLE) of the parameters \$\mu, \sigma, \pi\$, which are part of a gaussian mixture model. If you want to understand where the calculations come from, I think the easiest way is just to show you a short description of the algorithm and the formulas below.
Suppose the initial parameter vector is \$\hat\theta = (\{\hat\pi_l\}_{l=1}^L, \{\hat\mu_l\}_{l=1}^L, \{\hat\sigma_l^2\}_{l=1}^L)\$
Expectation step: Compute the expected value of the indicator variable that indicates which population (e.g., the population of skilled or unskilled managers) \$\alpha_{ij}\$ is drawn from:
\$\begin{align*}\hat{p}_{ijl} & = \hat{Pr}(\alpha_{ij}\text{ comes from Group }l) \\ & = \frac{\hat\pi_l\phi(\alpha_{ij};\hat\mu_l,\hat\sigma^2_l)}{\sum\limits_{l=1}^{L}\hat\pi_l\phi(\alpha_{ij};\hat\mu_l,\hat\sigma^2_l)}, i = 1, \ldots , N; j = 1, \ldots, M; l = 1, \ldots, L \end{align*}\$
where \$\phi(\cdot;\mu,\sigma^2)\$ is the density of the normal distribution \$\mathcal{N}(\mu, \sigma^2)\$Maximization Step: Compute the weighted means and variances with weights obtiained from the Expectation Step:
\$\tilde\mu_l=\frac{\sum_{ij}\hat{p}_{ijl}~\alpha_{ij}}{\sum_{ij}\hat{p}_{ijl}}, \tilde\sigma^2_l = \frac{\sum_{ij}\hat{p}_{ijl}~(\alpha_{ij} - \tilde\mu_l)^2}{\sum_{ij}\hat{p}_{ijl}}, \tilde\pi_l = \frac{\sum_{ij}\hat{p}_{ijl}}{MN}, l = 1, \ldots, L.\$Iterate between the Expectation Step and the Maximization Step until convergence.
