diff --git a/src/mdevaluate/correlation.py b/src/mdevaluate/correlation.py
index 7bcd3be..11be7da 100644
--- a/src/mdevaluate/correlation.py
+++ b/src/mdevaluate/correlation.py
@@ -199,7 +199,7 @@ def msd(
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
 
-def isf(
+def isf_raw(
     start_frame: CoordinateFrame,
     end_frame: CoordinateFrame,
     q: float = 22.7,
@@ -216,29 +216,59 @@ def isf(
         displacements = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
         distance = (displacements**2).sum(axis=1) ** 0.5
-        return np.sinc(distance * q / np.pi).mean()
+        return np.sinc(distance * q / np.pi)
     elif axis == "xy" or axis == "yx":
         distance = (displacements[:, [0, 1]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q)).mean()
+        return np.real(jn(0, distance * q))
     elif axis == "xz" or axis == "zx":
         distance = (displacements[:, [0, 2]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q)).mean()
+        return np.real(jn(0, distance * q))
     elif axis == "yz" or axis == "zy":
         distance = (displacements[:, [1, 2]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q)).mean()
+        return np.real(jn(0, distance * q))
     elif axis == "x":
         distance = np.abs(displacements[:, 0])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     elif axis == "y":
         distance = np.abs(displacements[:, 1])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     elif axis == "z":
         distance = np.abs(displacements[:, 2])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     else:
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
 
+def isf(
+    start_frame: CoordinateFrame,
+    end_frame: CoordinateFrame,
+    q: float = 22.7,
+    trajectory: Coordinates = None,
+    axis: str = "all",
+) -> float:
+    """
+    Incoherent intermediate scattering function averaged over all particles.
+    See isf_raw for details.
+    """
+    return isf_raw(start_frame, end_frame, q=q, trajectory=trajectory, axis=axis).mean()
+
+
+def isf_mean_var(
+    start_frame: CoordinateFrame,
+    end_frame: CoordinateFrame,
+    q: float = 22.7,
+    trajectory: Coordinates = None,
+    axis: str = "all",
+) -> float:
+    """
+    Incoherent intermediate scattering function averaged over all particles and the
+    variance.
+    See isf_raw for details.
+    """
+    values = isf_raw(start_frame, end_frame, q=q, trajectory=trajectory, axis=axis)
+    return values.mean(), values.var()
+
+
 def rotational_autocorrelation(
     start_frame: CoordinateFrame, end_frame: CoordinateFrame, order: int = 2
 ) -> float:
@@ -464,5 +494,8 @@ def non_gaussian_parameter(
         dimensions = 1
     else:
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
-
-    return (np.mean(r**2) / ((1 + 2 / dimensions) * (np.mean(r) ** 2))) - 1
+    
+    mean_r = np.mean(r)
+    if mean_r == 0.0:
+        return 0.0
+    return (np.mean(r**2) / ((1 + 2 / dimensions) * (mean_r ** 2))) - 1