code: def covariance_matrix(self, best_fit_params):
          
          cov_matrix = self._covariance_matrix(best_fit_params,
                                               objective=self.objective)
          if cov_matrix is None:
              # If the covariance matrix could not be computed we try again by
              # approximating the hessian with the jacobian.
  
              # VectorLeastSquares should be turned into a LeastSquares for
              # cov matrix calculation
              if self.objective.__class__ is VectorLeastSquares:
                  base = LeastSquares
              else:
                  base = self.objective.__class__
  
              class HessApproximation(base, HessianObjectiveJacApprox):
                  """
                  Class which impersonates ``base``, but which returns zeros
                  for the models Hessian. This will effectively result in the
                  calculation of the approximate Hessian by calculating
                  outer(J.T, J) when calling ``base.eval_hessian``.
                  """
  
              objective = HessApproximation(self.objective.model,
                                            self.objective.data)
              cov_matrix = self._covariance_matrix(best_fit_params,
                                                   objective=objective)
  
          return cov_matrix
docs: Given best fit parameters, this function finds the covariance matrix.
          This matrix gives the (co)variance in the parameters.
  
          :param best_fit_params: ``dict`` of best fit parameters as given by .best_fit_params()
          :return: covariance matrix.