def _gmm(profiles, n_groups, threshold, init='random', n_init=100):
    """Apply the classic Gaussian Mixture Model (GMM) to discover 
    resource groups [1]_.

    Parameters
    ----------
    profiles : DataFrame
        Constructed resource profiles.
    n_groups : int
        Expected number of resource groups.
    threshold : float
        A given threshold value in range [0, 1.0] for producing
        determined clustering from the fuzzy clustering results from GMM.
    init : {'random', 'mja', 'ahc', 'kmeans', 'plain'}, optional,         default 'random'
        Options for deciding the strategy for initialization. Could be 
        one of the following:

            - ``'random'``, use random initialization.
            - ``'mja'``, use the mining method Metric based on Joint 
              Activities as initialization method.
            - ``'ahc'``, use the hierarchical method of Agglomerative
              Hierarchical Clustering as initialization method.
            - ``'kmeans'``, use the classic clustering algorithm kMeans 
              as initialization method.
            - ``'plain'``, simply put resources into groups by in the 
              order of their ids.

        Note that if an option other than ``'random'`` is specified, then
        the initialization is performed once only.
    n_init : int, optional, default 100
        Number of times of random initialization (if specified) performs
        before training the clustering model.

    Returns
    -------
    list of frozensets
        Discovered resource groups.
    
    Raises
    ------
    ValueError
        If the specified option for initialization is invalid.
    
    See Also
    --------
    sklearn.mixture.GaussianMixture
    orgminer.OrganizationalModelMiner.community.graph_partitioning.mja
    orgminer.OrganizationalModelMiner.clustering.hierarchical.ahc
    sklearn.cluster.KMeans

    References
    ----------
    .. [1] Yang, J., Ouyang, C., Pan, M., Yu, Y., & ter Hofstede, A.
       (2018). Finding the Liberos: Discover Organizational Models with
       Overlaps. In *Proceedings of the 16th International Conference on
       Business Process Management*, pp. 339-355. Springer, Cham.
       `<https://doi.org/10.1007/978-3-319-98648-7_20>`_
    """
    print('Applying overlapping clustering -based GMM:')
    if init in ['mja', 'ahc', 'kmeans', 'plain']:
        warm_start = True
        if init == 'mja':
            from orgminer.OrganizationalModelMiner.community.graph_partitioning import _mja
            init_groups = _mja(profiles, n_groups)
        elif init == 'ahc':
            from .hierarchical import _ahc
            init_groups, _ = _ahc(profiles, n_groups)
        elif init == 'kmeans':
            init_groups = list(set() for i in range(n_groups))
            from sklearn.cluster import KMeans
            labels = KMeans(n_clusters=n_groups, random_state=0).fit_predict(
                profiles)
            for i, r in enumerate(sorted(profiles.index)):
                init_groups[labels[i]].add(r)
        else:
            init_groups = list(set() for i in range(n_groups))
            for i, r in enumerate(sorted(profiles.index)):
                init_groups[i % n_groups].add(r)
        print('Initialization done using {}:'.format(init))
    elif init == 'random':
        warm_start = False
    else:
        raise ValueError('Invalid value for parameter `{}`: {}'.format(
            'init', init))
    from sklearn.mixture import GaussianMixture
    if warm_start:
        from numpy import mean
        init_means = list()
        for g in init_groups:
            init_means.append(mean(profiles.loc[list(g)].values, axis=0))
        gmm_model = GaussianMixture(n_components=n_groups, covariance_type=
            'tied', tol=1e-09, max_iter=1000, n_init=1, random_state=0,
            means_init=init_means).fit(profiles.values)
    else:
        gmm_model = GaussianMixture(n_components=n_groups, covariance_type=
            'tied', tol=1e-09, max_iter=1000, n_init=n_init, init_params=
            'random').fit(profiles.values)
    posterior_pr = gmm_model.predict_proba(profiles.values)
    from numpy import nonzero, median, unique, count_nonzero, amin, percentile
    """
    threshold_ub = 0
    for pr in sorted(unique(posterior_pr), reverse=True):
        mbr_mat = posterior_pr >= pr
        if count_nonzero(mbr_mat.any(axis=1)) < len(profiles):
            pass
        else:
            if count_nonzero(mbr_mat.any(axis=0)) == n_groups:
                threshold_ub = pr
                break

    if threshold_ub != 0.0:
        threshold = amin([percentile(posterior_pr, 75), threshold_ub])
        membership_total = posterior_pr >= threshold
    else:
        membership_total = posterior_pr > 0.0
    """
    membership_total = posterior_pr >= threshold
    from collections import defaultdict
    groups = defaultdict(set)
    for i, membership in enumerate(membership_total):
        for j in nonzero(membership)[0]:
            groups[j].add(profiles.index[i])
    return [frozenset(g) for g in groups.values()]
