# imports
import numpy as np
from sklearn.datasets import load_iris
import pandas as pd
from sklearn.mixture import GaussianMixture as GM

clusters = [i for i in range(2,8)]
models = []
random_state = 42
n_init = 10

# train one model per cluster size
cols = ["petal width (cm)","petal length (cm)"]
models = [GM(c, n_init=n_init, init_params="random", random_state=random_state).\
          fit(X[cols]) for c in clusters]

# get best model from BIC
bic = [model.bic(X[cols]) for model in models]
metrics = pd.DataFrame({"clusters":clusters,"bic":bic})
metrics["best"] = metrics["bic"] == min(metrics["bic"])
metrics

bics = [model._n_parameters() * np.log(X.shape[0]) -2*model.score(X[cols]) for k,model in enumerate(models)]
bics