Bug fix in default computation of ranges in contour plots (#26)

Co-authored-by: Marina Evers <[email protected]>

uadapy/plotting/plots2D.py

@@ -107,7 +107,7 @@ def plot_contour(distributions,
     ranges : list or None, optional
         The ranges for the x and y axes. If None, the ranges are calculated based on the distributions.
     quantiles : list or None, optional
-        List of quantiles to use for determining isovalues. If None, the 99.7%, 95%, and 68% quantiles are used.
+        List of quantiles to use for determining isovalues. If None, the 95%, 75%, and 25% quantiles are used.
     seed : int
         Seed for the random number generator for reproducibility. It defaults to 55 if not provided.
     distrib_colors : list or None, optional
@@ -150,18 +150,36 @@ def plot_contour(distributions,
             distrib_colors.extend(additional_colors)
         palette = distrib_colors
 
+    # Determine default quantiles: 25%, 75%, and 95%
+    if quantiles is None:
+        quantiles = [25, 75, 95]
+    largest_quantile = max(quantiles)
+
     if ranges is None:
         min_val = np.zeros(distributions[0].mean().shape)+1000
         max_val = np.zeros(distributions[0].mean().shape)-1000
         cov_max = np.zeros(distributions[0].mean().shape)
+
+        # Dynamically set the expansion factor based on the largest quantile
+        if largest_quantile >= 99.99:
+            ef = 6  # 6 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99.9:
+            ef = 5  # 5 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99:
+            ef = 4  # 4 standard deviations for 99% quantiles
+        else:
+            ef = 3  # Default to 3 standard deviations
+
         for d in distributions:
             min_val=np.min(np.array([d.mean(), min_val]), axis=0)
             max_val=np.max(np.array([d.mean(), max_val]), axis=0)
             cov_max = np.max(np.array([np.diagonal(d.cov()), cov_max]), axis=0)
         cov_max = np.sqrt(cov_max)
-        ranges = [(mi-3*co, ma+3*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+        ranges = [(mi-ef*co, ma+ef*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+
     range_x = ranges[0]
     range_y = ranges[1]
+
     for i, d in enumerate(distributions):
         x = np.linspace(range_x[0], range_x[1], resolution)
         y = np.linspace(range_y[0], range_y[1], resolution)
@@ -178,18 +196,13 @@ def plot_contour(distributions,
         samples = d.sample(n_samples, seed)
         densities = d.pdf(samples)
         densities.sort()
-        if quantiles is None:
-            isovalues.append(densities[int((1 - 99.7/100) * n_samples)]) # 99.7% quantile
-            isovalues.append(densities[int((1 - 95/100) * n_samples)]) # 95% quantile
-            isovalues.append(densities[int((1 - 68/100) * n_samples)]) # 68% quantile
-        else:
-            quantiles.sort(reverse=True)
-            for quantile in quantiles:
-                if not 0 < quantile < 100:
-                    raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
-                elif int((1 - quantile/100) * n_samples) >= n_samples:
-                    raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
-                isovalues.append(densities[int((1 - quantile/100) * n_samples)])
+        quantiles.sort(reverse=True)
+        for quantile in quantiles:
+            if not 0 < quantile < 100:
+                raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
+            elif int((1 - quantile/100) * n_samples) >= n_samples:
+                raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
+            isovalues.append(densities[int((1 - quantile/100) * n_samples)])
 
         plt.contour(xv, yv, pdf, levels=isovalues, colors = [color])
 
@@ -224,7 +237,7 @@ def plot_contour_bands(distributions,
     ranges : list or None, optional
         The ranges for the x and y axes. If None, the ranges are calculated based on the distributions.
     quantiles : list or None, optional
-        List of quantiles to use for determining isovalues. If None, the 99.7%, 95%, and 68% quantiles are used.
+        List of quantiles to use for determining isovalues. If None, the 95%, 75%, and 25% quantiles are used.
     seed : int
         Seed for the random number generator for reproducibility. It defaults to 55 if not provided.
     show_plot : bool, optional
@@ -251,18 +264,36 @@ def plot_contour_bands(distributions,
     alpha_values = np.linspace(1/n_quantiles, 1.0, n_quantiles)  # Creates alpha values from 1/n to 1.0
     custom_cmap = utils.create_shaded_set2_colormap(alpha_values)
 
+    # Determine default quantiles: 25%, 75%, and 95%
+    if quantiles is None:
+        quantiles = [25, 75, 95]
+    largest_quantile = max(quantiles)
+
     if ranges is None:
         min_val = np.zeros(distributions[0].mean().shape)+1000
         max_val = np.zeros(distributions[0].mean().shape)-1000
         cov_max = np.zeros(distributions[0].mean().shape)
+
+        # Dynamically set the expansion factor based on the largest quantile
+        if largest_quantile >= 99.99:
+            ef = 6  # 6 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99.9:
+            ef = 5  # 5 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99:
+            ef = 4  # 4 standard deviations for 99% quantiles
+        else:
+            ef = 3  # Default to 3 standard deviations
+
         for d in distributions:
             min_val=np.min(np.array([d.mean(), min_val]), axis=0)
             max_val=np.max(np.array([d.mean(), max_val]), axis=0)
             cov_max = np.max(np.array([np.diagonal(d.cov()), cov_max]), axis=0)
         cov_max = np.sqrt(cov_max)
-        ranges = [(mi-3*co, ma+3*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+        ranges = [(mi-ef*co, ma+ef*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+
     range_x = ranges[0]
     range_y = ranges[1]
+
     for i, d in enumerate(distributions):
         x = np.linspace(range_x[0], range_x[1], resolution)
         y = np.linspace(range_y[0], range_y[1], resolution)
@@ -278,19 +309,14 @@ def plot_contour_bands(distributions,
         samples = d.sample(n_samples, seed)
         densities = d.pdf(samples)
         densities.sort()
-        if quantiles is None:
-            isovalues.append(densities[int((1 - 99.7/100) * n_samples)]) # 99.7% quantile
-            isovalues.append(densities[int((1 - 95/100) * n_samples)]) # 95% quantile
-            isovalues.append(densities[int((1 - 68/100) * n_samples)]) # 68% quantile
-        else:
-            quantiles.sort(reverse=True)
-            for quantile in quantiles:
-                if not 0 < quantile < 100:
-                    raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
-                elif int((1 - quantile/100) * n_samples) >= n_samples:
-                    raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
-                isovalues.append(densities[int((1 - quantile/100) * n_samples)])
-        isovalues.append(densities[-1])  # Minimum density value
+
+        quantiles.sort(reverse=True)
+        for quantile in quantiles:
+            if not 0 < quantile < 100:
+                raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
+            elif int((1 - quantile/100) * n_samples) >= n_samples:
+                raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
+            isovalues.append(densities[int((1 - quantile/100) * n_samples)])
 
         # Extract the subset of colors corresponding to the current Set2 color and its 3 alpha variations
         start_idx = i * n_quantiles

uadapy/plotting/plotsND.py

@@ -117,7 +117,7 @@ def plot_contour(distributions,
     ranges : list or None, optional
         Array of ranges for all dimensions. If None, the ranges are calculated based on the distributions.
     quantiles : list or None, optional
-        List of quantiles to use for determining isovalues. If None, the 99.7%, 95%, and 68% quantiles are used.
+        List of quantiles to use for determining isovalues. If None, the 95%, 75%, and 25% quantiles are used.
     seed : int
         Seed for the random number generator for reproducibility. It defaults to 55 if not provided.
     distrib_colors : list or None, optional
@@ -147,6 +147,11 @@ def plot_contour(distributions,
     if isinstance(distributions, Distribution):
         distributions = [distributions]
 
+    # Determine default quantiles: 25%, 75%, and 95%
+    if quantiles is None:
+        quantiles = [25, 75, 95]
+    largest_quantile = max(quantiles)
+
     # Generate colors
     if distrib_colors is None:
         if colorblind_safe:
@@ -168,12 +173,24 @@ def plot_contour(distributions,
         min_val = np.zeros(distributions[0].mean().shape)+1000
         max_val = np.zeros(distributions[0].mean().shape)-1000
         cov_max = np.zeros(np.diagonal(distributions[0].cov()).shape)
+
+        # Dynamically set the expansion factor based on the largest quantile
+        if largest_quantile >= 99.99:
+            ef = 6  # 6 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99.9:
+            ef = 5  # 5 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99:
+            ef = 4  # 4 standard deviations for 99% quantiles
+        else:
+            ef = 3  # Default to 3 standard deviations
+
         for d in distributions:
             min_val=np.min(np.array([d.mean(), min_val]), axis=0)
             max_val=np.max(np.array([d.mean(), max_val]), axis=0)
             cov_max = np.max(np.array([np.diagonal(d.cov()), cov_max]), axis=0)
         cov_max = np.sqrt(cov_max)
-        ranges = [(mi-3*co, ma+3*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+        ranges = [(mi-ef*co, ma+ef*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+
     fig, axes = plt.subplots(nrows=n_dims, ncols=n_dims)
     for i, ax in enumerate(axes.flat):
         # Hide all ticks and labels
@@ -200,18 +217,13 @@ def plot_contour(distributions,
         samples = d.sample(n_samples, seed)
         densities = d.pdf(samples)
         densities.sort()
-        if quantiles is None:
-            isovalues.append(densities[int((1 - 99.7/100) * n_samples)]) # 99.7% quantile
-            isovalues.append(densities[int((1 - 95/100) * n_samples)]) # 95% quantile
-            isovalues.append(densities[int((1 - 68/100) * n_samples)]) # 68% quantile
-        else:
-            quantiles.sort(reverse=True)
-            for quantile in quantiles:
-                if not 0 < quantile < 100:
-                    raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
-                elif int((1 - quantile/100) * n_samples) >= n_samples:
-                    raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
-                isovalues.append(densities[int((1 - quantile/100) * n_samples)])
+        quantiles.sort(reverse=True)
+        for quantile in quantiles:
+            if not 0 < quantile < 100:
+                raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
+            elif int((1 - quantile/100) * n_samples) >= n_samples:
+                raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
+            isovalues.append(densities[int((1 - quantile/100) * n_samples)])
 
         for i, j in zip(*np.triu_indices_from(axes, k=1)):
             for x, y in [(i, j), (j, i)]:
@@ -271,7 +283,7 @@ def plot_contour_samples(distributions,
     ranges : list or None, optional
         Array of ranges for all dimensions. If None, the ranges are calculated based on the distributions.
     quantiles : list or None, optional
-        List of quantiles to use for determining isovalues. If None, the 99.7%, 95%, and 68% quantiles are used.
+        List of quantiles to use for determining isovalues. If None, the 95%, 75%, and 25% quantiles are used.
     seed : int
         Seed for the random number generator for reproducibility. It defaults to 55 if not provided.
     distrib_colors : list or None, optional
@@ -301,6 +313,11 @@ def plot_contour_samples(distributions,
     if isinstance(distributions, Distribution):
         distributions = [distributions]
 
+    # Determine default quantiles: 25%, 75%, and 95%
+    if quantiles is None:
+        quantiles = [25, 75, 95]
+    largest_quantile = max(quantiles)
+
     # Generate colors
     if distrib_colors is None:
         if colorblind_safe:
@@ -322,12 +339,24 @@ def plot_contour_samples(distributions,
         min_val = np.zeros(distributions[0].mean().shape)+1000
         max_val = np.zeros(distributions[0].mean().shape)-1000
         cov_max = np.zeros(np.diagonal(distributions[0].cov()).shape)
+
+        # Dynamically set the expansion factor based on the largest quantile
+        if largest_quantile >= 99.99:
+            ef = 6  # 6 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99.9:
+            ef = 5  # 5 standard deviations for 99.9% quantiles
+        elif largest_quantile >= 99:
+            ef = 4  # 4 standard deviations for 99% quantiles
+        else:
+            ef = 3  # Default to 3 standard deviations
+
         for d in distributions:
             min_val=np.min(np.array([d.mean(), min_val]), axis=0)
             max_val=np.max(np.array([d.mean(), max_val]), axis=0)
             cov_max = np.max(np.array([np.diagonal(d.cov()), cov_max]), axis=0)
         cov_max = np.sqrt(cov_max)
-        ranges = [(mi-3*co, ma+3*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+        ranges = [(mi-ef*co, ma+ef*co) for mi,ma, co in zip(min_val, max_val, cov_max)]
+
     fig, axes = plt.subplots(nrows=n_dims, ncols=n_dims)
     for i, ax in enumerate(axes.flat):
         # Hide all ticks and labels
@@ -353,18 +382,13 @@ def plot_contour_samples(distributions,
         samples = d.sample(n_samples, seed)
         densities = d.pdf(samples)
         densities.sort()
-        if quantiles is None:
-            isovalues.append(densities[int((1 - 99.7/100) * n_samples)]) # 99.7% quantile
-            isovalues.append(densities[int((1 - 95/100) * n_samples)]) # 95% quantile
-            isovalues.append(densities[int((1 - 68/100) * n_samples)]) # 68% quantile
-        else:
-            quantiles.sort(reverse=True)
-            for quantile in quantiles:
-                if not 0 < quantile < 100:
-                    raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
-                elif int((1 - quantile/100) * n_samples) >= n_samples:
-                    raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
-                isovalues.append(densities[int((1 - quantile/100) * n_samples)])
+        quantiles.sort(reverse=True)
+        for quantile in quantiles:
+            if not 0 < quantile < 100:
+                raise ValueError(f"Invalid quantile: {quantile}. Quantiles must be between 0 and 100 (exclusive).")
+            elif int((1 - quantile/100) * n_samples) >= n_samples:
+                raise ValueError(f"Quantile {quantile} results in an index that is out of bounds.")
+            isovalues.append(densities[int((1 - quantile/100) * n_samples)])
 
         for i, j in zip(*np.triu_indices_from(axes, k=1)):
             for x, y in [(i, j), (j, i)]: