code: def _likelihood(self, chol_sigma_s_rhos, log_det_psi, chol_sigma_s,
                      trace_xt_invsigma2_x, inv_sigma_s_rhos, wt_invpsi_x,
                      samples):
          
          log_det = (np.log(np.diag(chol_sigma_s_rhos) ** 2).sum() + log_det_psi
                     + np.log(np.diag(chol_sigma_s) ** 2).sum())
          loglikehood = -0.5 * samples * log_det - 0.5 * trace_xt_invsigma2_x
          loglikehood += 0.5 * np.trace(
              wt_invpsi_x.T.dot(inv_sigma_s_rhos).dot(wt_invpsi_x))
          # + const --> -0.5*nTR*nvoxel*subjects*math.log(2*math.pi)
  
          return loglikehood
docs: Calculate the log-likelihood function
  
  
          Parameters
          ----------
  
          chol_sigma_s_rhos : array, shape=[features, features]
              Cholesky factorization of the matrix (Sigma_S + sum_i(1/rho_i^2)
              * I)
  
          log_det_psi : float
              Determinant of diagonal matrix Psi (containing the rho_i^2 value
              voxels_i times).
  
          chol_sigma_s : array, shape=[features, features]
              Cholesky factorization of the matrix Sigma_S
  
          trace_xt_invsigma2_x : float
              Trace of :math:`\\sum_i (||X_i||_F^2/\\rho_i^2)`
  
          inv_sigma_s_rhos : array, shape=[features, features]
              Inverse of :math:`(\\Sigma_S + \\sum_i(1/\\rho_i^2) * I)`
  
          wt_invpsi_x : array, shape=[features, samples]
  
          samples : int
              The total number of samples in the data.
  
  
          Returns
          -------
  
          loglikehood : float
              The log-likelihood value.