diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index ea478462e860f265d812d7627bb78b689fbe0a7e..becf16977da432623ef27b487782ba3c40414068 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -4,10 +4,10 @@
 
 TODO: Give option to choose from multiwavelet degrees (first, last or
     highest magnitude) -> Done
-TODO: Change method of calculating quartiles
-TODO: Include overlapping cells in quartile calculation (if needed)
+TODO: Change method of calculating quartiles -> Done
+TODO: Include overlapping cells in quartile calculation (if needed) -> Done
 TODO: Determine max_value for Theoretical only over highest degree
-    -> Done (optional)
+    -> Done (now optional)
 TODO: Check if indexing in wavelets is correct
 TODO: Combine get_cells() and _get_cells()
 TODO: Add TC condition to only flag cell if left-adjacent one is flagged as
@@ -380,7 +380,9 @@ class Boxplot(WaveletDetector):
         Flag whether outer fences should be adjusted using global mean.
     extreme_outlier_only : bool
         Flag whether outliers also have to be detected in neighbouring folds.
-    folds : ndarray
+    quantile_method : str
+        Method used to calculate quantiles.
+    fold_indices : ndarray
         Array with indices for elements of each fold (including
         overlaps).
 
@@ -405,6 +407,7 @@ class Boxplot(WaveletDetector):
         if self._mesh.num_grid_cells < self._fold_len:
             self._fold_len = self._mesh.num_grid_cells
 
+        self._quantile_method = config.pop('quantile_method', 'weibull')
         num_overlapping_cells = config.pop('num_overlapping_cells', 1)
         num_folds = self._mesh.num_grid_cells//self._fold_len
         self._fold_indices = np.zeros([num_folds,
@@ -432,21 +435,12 @@ class Boxplot(WaveletDetector):
             List of indices of all detected troubled cells.
 
         """
-        # Select and sort fold domains
+        # Determine quartiles of folds
         folds = coeffs[self._fold_indices]
-        folds.sort()
-
-        # Determine quartile parameters
-        boundary_index = self._fold_len//4
-        balance_factor = self._fold_len/4.0 - boundary_index
-
-        # Determine first and third quartiles
-        first_quartiles = (1-balance_factor) \
-            * folds[:, boundary_index-1] \
-            + balance_factor * folds[:, boundary_index]
-        third_quartiles = (1-balance_factor) \
-            * folds[:, 3*boundary_index-1]\
-            + balance_factor * folds[:, 3*boundary_index]
+        first_quartiles = np.quantile(folds, 0.25, axis=1,
+                                      method=self._quantile_method)
+        third_quartiles = np.quantile(folds, 0.75, axis=1,
+                                      method=self._quantile_method)
 
         # Determine bounds based on quartiles of a boxplot
         lower_bounds = np.zeros(len(first_quartiles) + 2)