Gaussian Model

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal as multiNorm

def draw_vector(p,v):
    # Plots the vector v with  origin at point p
    plt.quiver(p[0],p[1], v[0], v[1])
    plt.axis('equal')

# Datos Iniciales.
X = list(map(np.array, ((-4,0), (2,0), (-2,6)))) # Puntos como columnas.
Gaussians = [[0.8, np.array((-1,0)), np.matrix(((1,0), (0,1)))], [0.2, np.array((4,0)), np.matrix(((2,-1), (-1,2)))]]
K = len(Gaussians)

# Expectation Step.
R = [[], []] # [R1, R2]
for k in range(K):
    for x in X:
        Numerator = Gaussians[k][0] * multiNorm.pdf(x, Gaussians[k][1], Gaussians[k][2])
        Denominator = sum([Gaussians[l][0]*multiNorm.pdf(x, Gaussians[l][1], Gaussians[l][2]) for l in range(K)])
        R[k].append(Numerator/Denominator)

print('Responsabilities:')
for k in range(K):
    for n in range(len(X)):
        print(f'K{k+1}N{n+1}', R[k][n])

# Maximisation Step.
Xc = [x.T for x in list(map(np.matrix, ((-4,0), (2,0), (-2,6))))] # Puntos como columnas.
UpdatedGaussians = [[], []]

for k in range(K):
    N = sum([R for R in R[k]])
    Pi = N/len(Xc)
    Mean = (1/N) * sum([R[k][n]*Xc[n] for n in range(len(Xc))])
    Cov = (1/N) * sum([R[k][n] * (Xc[n]-Mean) * (Xc[n]-Mean).T for n in range(len(Xc))])
    UpdatedGaussians[k] = [Pi, Mean, Cov]

print('New Gaussians Parameters:', '\n')
for k in range(K):
    print(f'π{k+1}',UpdatedGaussians[k][0])
    print(f'μ{k+1}',UpdatedGaussians[k][1])
    print(f'Σ{k+1}',UpdatedGaussians[k][2])
    print('\n')

# Convergence.
n = sum([np.log(sum([Gaussians[k][0] * multiNorm.pdf(x, Gaussians[k][1], Gaussians[k][2]) for k in range(K)])) for x in X])
print(f'Log-Likehood Value:{n}')

# Draw Models.
for k in range(K):
    Center = UpdatedGaussians[k][1]
    Directions, Eigenvalues = np.linalg.svd(UpdatedGaussians[k][2])[0], np.linalg.svd(UpdatedGaussians[k][2])[1]
    for i in range (len(Directions)): draw_vector(Center, np.squeeze(np.asarray(Directions[:,i])))