code: def three_classes_GMM_fit(
      x,
      test=None,
      alpha=0.01,
      prior_strength=100,
      verbose=0,
      fixed_scale=False,
      mpaxes=None,
      bias=0,
      theta=0,
      return_estimator=False,
  ):
      """Fit the data with a 3-classes Gaussian Mixture Model,
      i.e. compute some probability that the voxels of a certain map
      are in class disactivated, null or active
  
      Parameters
      ----------
      x: array of shape (nvox,1)
        The map to be analysed
      test: array of shape(nbitems,1), optional
        the test values for which the p-value needs to be computed
        by default (if None), test=x
      alpha: float, optional
        the prior weights of the positive and negative classes
      prior_strength: float, optional
        the confidence on the prior (should be compared to size(x))
      verbose: int
        verbosity mode
      fixed_scale: bool, optional
        boolean, variance parameterization. if True, the variance is locked to 1
        otherwise, it is estimated from the data
      mpaxes:
        axes handle used to plot the figure in verbose mode
        if None, new axes are created
      bias:  bool
        allows a rescaling of the posterior probability
        that takes into account the threshold theta. Not rigorous.
      theta: float
        the threshold used to correct the posterior p-values
        when bias=1; normally, it is such that test>theta
        note that if theta = -np.inf, the method has a standard behaviour
      return_estimator: boolean, optional
        If return_estimator is true, the estimator object is
        returned.
  
      Returns
      -------
      bfp : array of shape (nbitems,3):
          the posterior probability of each test item belonging to each component
          in the GMM (sum to 1 across the 3 classes)
          if np.size(test)==0, i.e. nbitem==0, None is returned
      estimator : nipy.labs.clustering.GMM object
          The estimator object, returned only if return_estimator is true.
  
      Notes
      -----
      Our convention is that:
  
      * class 1 represents the negative class
      * class 2 represents the null class
      * class 3 represents the positive class
      """
      from ..clustering.bgmm import VBGMM
      from ..clustering.gmm import GridDescriptor
  
      nvox = np.size(x)
      x = np.reshape(x, (nvox, 1))
      if test is None:
          test = x
      if np.size(test) == 0:
          return None
  
      sx = np.sort(x, 0)
      nclasses = 3
  
      # set the priors from a reasonable model of the data (!)
      # prior means
      mb0 = np.mean(sx[: int(alpha * nvox)])
      mb2 = np.mean(sx[int((1 - alpha) * nvox) :])
      prior_means = np.reshape(np.array([mb0, 0, mb2]), (nclasses, 1))
      if fixed_scale:
          prior_scale = (
              np.ones((nclasses, 1, 1)) * 1.0 / (prior_strength)
          )
      else:
          prior_scale = (
              np.ones((nclasses, 1, 1))
              * 1.0
              / (prior_strength * np.var(x))
          )
      prior_dof = np.ones(nclasses) * prior_strength
      prior_weights = (
          np.array([alpha, 1 - 2 * alpha, alpha]) * prior_strength
      )
      prior_shrinkage = np.ones(nclasses) * prior_strength
  
      # instantiate the class and set the priors
      BayesianGMM = VBGMM(
          nclasses,
          1,
          prior_means,
          prior_scale,
          prior_weights,
          prior_shrinkage,
          prior_dof,
      )
      BayesianGMM.set_priors(
          prior_means,
          prior_weights,
          prior_scale,
          prior_dof,
          prior_shrinkage,
      )
  
      # estimate the model
      BayesianGMM.estimate(x, delta=1.0e-8, verbose=max(0, verbose - 1))
  
      # create a sampling grid
      if verbose or bias:
          gd = GridDescriptor(1)
          gd.set([x.min(), x.max()], 100)
          gdm = gd.make_grid().squeeze()
          lj = BayesianGMM.likelihood(gd.make_grid())
  
      # estimate the prior weights
      bfp = BayesianGMM.likelihood(test)
      if bias:
          lw = np.sum(lj[gdm > theta], 0)
          weights = BayesianGMM.weights / (BayesianGMM.weights.sum())
          bfp = (lw / weights) * BayesianGMM.slikelihood(test)
  
      if verbose and (mpaxes is not False):
          BayesianGMM.show_components(x, gd, lj, mpaxes)
  
      bfp = (bfp.T / bfp.sum(1)).T
      if not return_estimator:
          return bfp
      else:
          return bfp, BayesianGMM
  
file_path: nipy/nipy/algorithms/statistics/empirical_pvalue.py
imports: import numpy
  import scipy.stats
  from numpy.linalg import pinv
label_desc: Use image node_id 6 for calling a global function with example usage: three_classes_GMM_fit(x, test, alpha, prior_strength, verbose, fixed_scale, mpaxes, bias, theta, return_estimator) and returns: None, bfp, bfp,  BayesianGMM