diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index 69f919f7..6ec89a04 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -353,6 +353,11 @@ def func(a, x):
     return a[0] * x1 ** 2 + a[1] * x2
 ```
 
+`pyerrors` also supports correlated fits which can be triggered via the parameter `correlated_fit=True`.
+Details about how the required covariance matrix is estimated can be found in `pyerrors.obs.covariance`.
+
+Direct visualizations of the performed fits can be triggered via `resplot=True` or `qqplot=True`. For all available options see `pyerrors.fits.least_squares`.
+
 ## Total least squares fits
 `pyerrors` can also fit data with errors on both the dependent and independent variables using the total least squares method also referred to orthogonal distance regression as implemented in [scipy](https://docs.scipy.org/doc/scipy/reference/odr.html), see `pyerrors.fits.least_squares`. The syntax is identical to the standard least squares case, the only diffrence being that `x` also has to be a `list` or `numpy.array` of `Obs`.
 
diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index 6a428e65..c1925560 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -86,7 +86,7 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
         ```
 
         It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
-        will not work
+        will not work.
     priors : list, optional
         priors has to be a list with an entry for every parameter in the fit. The entries can either be
         Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like
@@ -95,24 +95,25 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs):
         If true all output to the console is omitted (default False).
     initial_guess : list
         can provide an initial guess for the input parameters. Relevant for
-                     non-linear fits with many parameters.
+        non-linear fits with many parameters.
     method : str, optional
         can be used to choose an alternative method for the minimization of chisquare.
         The possible methods are the ones which can be used for scipy.optimize.minimize and
         migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
         Reliable alternatives are migrad, Powell and Nelder-Mead.
-    resplot : bool
-        If true, a plot which displays fit, data and residuals is generated (default False).
-    qqplot : bool
-        If true, a quantile-quantile plot of the fit result is generated (default False).
-    expected_chisquare : bool
-        If true prints the expected chisquare which is
-        corrected by effects caused by correlated input data.
-        This can take a while as the full correlation matrix
-        has to be calculated (default False).
     correlated_fit : bool
-        If true, use the full correlation matrix in the definition of the chisquare
-        (only works for prior==None and when no method is given, at the moment).
+        If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
+        For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`.
+        In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix).
+        This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning).
+        At the moment this option only works for `prior==None` and when no `method` is given.
+    expected_chisquare : bool
+        If True estimates the expected chisquare which is
+        corrected by effects caused by correlated input data (default False).
+    resplot : bool
+        If True, a plot which displays fit, data and residuals is generated (default False).
+    qqplot : bool
+        If True, a quantile-quantile plot of the fit result is generated (default False).
     '''
     if priors is not None:
         return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 6889fd3f..8e4351a8 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -1348,7 +1348,9 @@ def covariance(obs, visualize=False, correlation=False, **kwargs):
 
     Notes
     -----
-    The covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. For observables defined on a single ensemble this is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
+    The covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
+    $$v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v_i\in\mathbb{R}^N$, while such an identity does not hold for larger windows/lags.
+    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
     $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
     This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
     '''