docs: def inference(self, in_dims, out_dims, value=None):
          """ Perform Bayesian inference on the gmm. Let's call V = V1...Vd the d-dimensional space on which the current GMM is defined, such that it represents P(V). Let's call X and Y to disjoint subspaces of V, with corresponding dimension indices in ran. This method returns the GMM for P(Y | X=value).
  
          :param list in_dims: the dimension indices of X (a subset of range(d)). This can be the empty list if one want to compute the marginal P(Y).
  
          :param list out_dims: the dimension indices of Y (a subset of range(d), without intersection with in_dims).
  
          :param numpy.array value: the value of X for which one want to compute the conditional (ignored of in_dims=[]).
  
          :returns: the gmm corresponding to P(Y | X=value) (or to P(Y) if in_dims=[])
  
          .. note:: For example, if X = V1...Vm and Y = Vm+1...Vd, then P(Y | X=v1...vm) is returned by self.inference(in_dims=range(m), out_dims=range(m, d), array([v1, ..., vm])).
          """
  
          if self.covariance_type != 'diag' and self.covariance_type != 'full':
              raise ValueError("covariance type other than 'full' and 'diag' not allowed")
          in_dims = array(in_dims)
          out_dims = array(out_dims)
          value = array(value)
  
          means = zeros((self.n_components, len(out_dims)))
          covars = zeros((self.n_components, len(out_dims), len(out_dims)))
          weights = zeros((self.n_components,))
          if in_dims.size:
              for k, (weight_k, mean_k, covar_k) in enumerate(self):
  
                  sig_in = covar_k[ix_(in_dims, in_dims)]
                  inin_inv = matrix(sig_in).I
                  out_in = covar_k[ix_(out_dims, in_dims)]
                  mu_in = mean_k[in_dims].reshape(-1, 1)
                  means[k, :] = (mean_k[out_dims] +
                                 (out_in *
                                  inin_inv *
                                  (value.reshape(-1, 1) - mu_in)).T)
                  if self.covariance_type == 'full':
                      covars[k, :, :] = (covar_k[ix_(out_dims, out_dims)] -
                                         out_in *
                                         inin_inv *
                                         covar_k[ix_(in_dims, out_dims)])
                  elif self.covariance_type == 'diag':
                      covars[k, :] = covar_k[out_dims]
  
                  weights[k] = weight_k * Gaussian(mu_in.reshape(-1,),
                                                   sig_in).normal(value.reshape(-1,))
              weights /= sum(weights)
          else:
              means = self.means_[:, out_dims]
              if self.covariance_type == 'full':
                  covars = self.covariances_[ix_(range(self.n_components), out_dims, out_dims)]
              if self.covariance_type == 'diag':
                  covars = self.covariances_[ix_(range(self.n_components), out_dims)]
              weights = self.weights_
  
          res = GMM(n_components=self.n_components,
                    covariance_type=self.covariance_type)
          res.weights_ = weights
          res.means_ = means
          res.covariances_ = covars
          return res