docs: def eval_hessian(self, ordered_parameters=[], **parameters):
          """
          Hessian for log-likelihood is defined as
          :math:`\\nabla^2_{\\vec{p}}( \\log( L(\\vec{p} | \\vec{x})))`.
  
          :param parameters: values for the fit parameters.
          :return: array of length number of ``Parameter``'s in the model, with all partial derivatives evaluated at p, data.
          """
          evaluated_func = super(LogLikelihood, self).__call__(
              ordered_parameters, **parameters
          )
          evaluated_jac = super(LogLikelihood, self).eval_jacobian(
              ordered_parameters, **parameters
          )
          evaluated_hess = super(LogLikelihood, self).eval_hessian(
              ordered_parameters, **parameters
          )
  
          result = 0
          for f, jac_comp, hess_comp in zip(evaluated_func, evaluated_jac, evaluated_hess):
              # Outer product
              jac_outer_jac = np.einsum('i...,j...->ij...', jac_comp, jac_comp)
              dd_logf = - hess_comp / f[np.newaxis, np.newaxis, ...] + \
                        (1 / f**2)[np.newaxis, np.newaxis, ...] * jac_outer_jac
              # We sum away everything except the matrices in the axes 0 & 1.
              axes = tuple(range(2, len(dd_logf.shape)))
              result += np.sum(dd_logf, axis=axes, keepdims=False)
          else:
              return np.atleast_2d(np.squeeze(np.array(result)))