model selection for GaussianMixture by using GridSearch
<p>I'd like to use the function GaussianMixture by scikit-learn, and I have to perform model selection.
I want to do it by using GridSearchCV, and I would like to use for the selection the BIC and the AIC.
Both these values are implemented into GaussianMixture(), but I don't know how to insert them into the definition of my custom scorer, since the function </p>

<pre><code>make_scorer(score_func, greater_is_better=True, needs_proba=False, needs_threshold=False, **kwargs)
</code></pre>

<p>that I am using to create my custom scorer takes as input a function score_funct, that has to be defined as </p>

<pre><code>score_func(y, y_pred, **kwargs)
</code></pre>

<p>Can someone help me?</p>

<p>Using the BIC/AIC is an <em>alternative</em> to using cross validation. <code>GridSearchCV</code> selects models using cross validation. To perform model selection using the BIC/AIC we have to do something a little different. Let's take an example where we generate samples from two Gaussians, and then try to fit them using scikit-learn. </p>

<pre><code>import numpy as np
X1 = np.random.multivariate_normal([0.,0.],[[1.,0.],[0.,1.]],10000)
X2 = np.random.multivariate_normal([10.,10.],[[1.,0.],[0.,1.]],10000)
X = np.vstack((X1,X2))
np.random.shuffle(X)
</code></pre>

<p><a href="http://i.stack.imgur.com/3SxTu.png" rel="nofollow"><img src="http://i.stack.imgur.com/3SxTu.png" alt="enter image description here"></a></p>

<p><strong>Method 1: Cross-validation</strong></p>

<p><a href="https://en.wikipedia.org/wiki/Cross-validation_%28statistics%29" rel="nofollow">Cross validation</a> involves splitting the data into pieces. One then fits the model on some of the pieces ('training') and tests how well it performs on the remaining pieces ('validating'). This guards against over-fitting. Here we will use two-fold cross validation, where we split the data in half.</p>

<pre><code>from sklearn.mixture import GaussianMixture
from sklearn.model_selection import GridSearchCV
import matplotlib.pyplot as plt

#check 1-&gt;4 components
tuned_parameters = {'n_components': np.array([1,2,3,4])}
#construct grid search object that uses 2 fold cross validation
clf = GridSearchCV(GaussianMixture(),tuned_parameters,cv=2)
#fit the data
clf.fit(X)
#plot the number of Gaussians against their rank
plt.scatter(clf.cv_results_['param_n_components'],\
            clf.cv_results_['rank_test_score'])
</code></pre>

<p>We can see that 2-fold cross validation favours two Gaussian components, as we expect.</p>

<p><a href="http://i.stack.imgur.com/oQ9cg.png" rel="nofollow"><img src="http://i.stack.imgur.com/oQ9cg.png" alt="enter image description here"></a>   </p>

<p><strong>Method 2: BIC/AIC</strong></p>

<p>Instead of using cross-validation, we can evaluate the <a href="https://en.wikipedia.org/wiki/Bayesian_information_criterion" rel="nofollow">BIC</a> using the best-fit model given each number of Gaussians. We then choose the model that has the lowest BIC. The procedure would be identical if one used the AIC (although it is a different statistic, and can provide different answers: but your code structure would be identical to below).</p>

<pre><code>bic = np.zeros(4)
n = np.arange(1,5)
models = []
#loop through each number of Gaussians and compute the BIC, and save the model
for i,j in enumerate(n):
    #create mixture model with j components
    gmm = GaussianMixture(n_components=j)
    #fit it to the data
    gmm.fit(X)
    #compute the BIC for this model
    bic[i] = gmm.bic(X)
    #add the best-fit model with j components to the list of models
    models.append(gmm)
</code></pre>

<p>After carrying out this procedure, we can plot the number of Gaussians against the BIC.</p>

<pre><code>plt.plot(n,bic)
</code></pre>

<p><a href="http://i.stack.imgur.com/SwZ2R.png" rel="nofollow"><img src="http://i.stack.imgur.com/SwZ2R.png" alt="enter image description here"></a></p>

<p>So we can see that the BIC is minimised for two Gaussians, so the best model 
according to this method also has two components.</p>

<p>Because I took 10000 samples from two very well-separated Gaussians (i.e. the distance between their centres is much larger than either of their dispersions), the answer was very clear-cut. This is not always the case, and often neither of these methods will confidently tell you which number of Gaussians to use, but rather some sensible range.</p>

python scikit-learn grid-search