diff --git a/quickstats/_version.py b/quickstats/_version.py
index 54ed38f7caa92621e4aaa2b85cad5b2d58c20f62..83c2f44e7f4117744a7ae074c109510edab08e5e 100644
--- a/quickstats/_version.py
+++ b/quickstats/_version.py
@@ -1 +1 @@
-__version__ = "0.8.3.5.6"
+__version__ = "0.8.3.5.7"
diff --git a/quickstats/plots/abstract_plot.py b/quickstats/plots/abstract_plot.py
index 854d7f1fc151dd9c63f4e7345a0da12d37fb9dc3..d66d136f7c354c441a0a6a8235da09d5979b7449 100644
--- a/quickstats/plots/abstract_plot.py
+++ b/quickstats/plots/abstract_plot.py
@@ -58,7 +58,8 @@ from .artists import (
     FillBetween,
     ErrorBand,
     Annotation,
-    Text
+    Text,
+    Ellipse
 )
 
 # Type variables for better type hints
@@ -441,6 +442,33 @@ class AbstractPlot(AbstractObject):
             raise ValueError(f'config option not set: {key}')
         return config[key]
 
+    def get_target_color(
+        self,
+        name: str,
+        target: Optional[str] = None,
+        fallback: bool = True
+    ) -> Optional[str]:
+        """
+        Get color for a domain.
+
+        Parameters
+        ----------
+        name : str
+            The name to get the color for
+        ...
+        fallback : bool, default False
+            Whether to fall back to name-only lookup
+
+        Returns
+        -------
+        Optional[str]
+            The target color if found
+        """
+        domain = self.color_map.format(target, name)
+        if domain not in self.color_map and fallback:
+            return self.color_map.get(name)
+        return self.color_map.get(domain)
+
     def get_target_label(
         self,
         name: str,
@@ -524,6 +552,19 @@ class AbstractPlot(AbstractObject):
         artist = FillBetween(x=x, y1=y1, y2=y2, label=label, styles=kwargs)
         self.add_artist(artist, name=name or label)
 
+    def add_ellipse(
+        self,
+        xy: Tuple[float, float],
+        width: float,
+        height: float,
+        angle: float = 0,
+        label: Optional[str] = None,
+        name: Optional[str] = None,
+        **kwargs
+    ) -> None:
+        artist = Ellipse(xy=xy, width=width, height=height, angle=angle, label=label, styles=kwargs)
+        self.add_artist(artist, name=name or label)
+
     def add_errorband(
         self,
         x: ArrayLike,
@@ -772,7 +813,7 @@ class AbstractPlot(AbstractObject):
         handles: List[Artist] = []
         labels: List[str] = []
 
-        targets = targets or self.legend_order
+        targets = targets or self.legend_order or self.get_labelled_legend_domains()
 
         try:
             for name in targets:
diff --git a/quickstats/plots/artists.py b/quickstats/plots/artists.py
index f263d6b61ca29ed007d46a271bc12adfe2120c92..80451e11fc7e744290419cfa64ffc2729dd72032 100644
--- a/quickstats/plots/artists.py
+++ b/quickstats/plots/artists.py
@@ -1,5 +1,7 @@
 from typing import Optional, Union, Any, Dict, Tuple, Callable
 
+from matplotlib.patches import Ellipse as Ellipse_mpl
+
 from quickstats.core import mappings as mp
 from quickstats.core.decorators import dataclass_ex
 from quickstats.core.typing import ArrayLike
@@ -136,6 +138,30 @@ class FillBetween(SingleArtist):
         )
         return handle
 
+@dataclass_ex(kw_only=True)
+class Ellipse(SingleArtist):
+    
+    xy: Tuple[float, float]
+    width: float
+    height: float
+    angle: float
+
+    def draw(
+        self,
+        ax,
+        base_styles_map: Optional[StylesMapType] = None
+    ):
+        styles = self.resolve_styles(base_styles_map, keys='ellipse')
+        handle = Ellipse_mpl(
+            xy=self.xy,
+            width=self.width,
+            height=self.height,
+            angle=self.angle,
+            **styles
+        )
+        ax.add_patch(handle)
+        return handle
+
 @dataclass_ex(kw_only=True)
 class ErrorBand(MixedArtist):
 
diff --git a/quickstats/plots/likelihood_2D_plot.py b/quickstats/plots/likelihood_2D_plot.py
index 00068138fe4d6cf661a7b24aa291a983f72aeb3b..88199342d5b42a8f479fb8e3c1356480e1fb13a4 100644
--- a/quickstats/plots/likelihood_2D_plot.py
+++ b/quickstats/plots/likelihood_2D_plot.py
@@ -34,14 +34,15 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
         'bestfit': {
             'marker': '*',
             'linewidth': 0,
-            'markersize': 15
+            'markersize': 20,
+            'markeredgecolor': 'black'
         },
         'point': {
             'linewidth': 0,
             'marker': '*',
             'markersize': 20,
             'color': '#E9F1DF',
-            'markeredgecolor': 'black'            
+            'markeredgecolor': 'black'
         },
         'contourf': {
             'extend': 'min'
@@ -77,7 +78,8 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
         },
         'remove_nan_points_within_distance': None,
         'shade_nan_points': False,
-        'alphashape_alpha': 0.1
+        'alphashape_alpha': 0.1,
+        'distinct_colors': False
     }
 
     def __init__(
@@ -128,7 +130,7 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
         ]:
             if domain in self.color_map:
                 color = self.color_map[domain]
-                if color not in used_colors:
+                if (color not in used_colors) or (not self.config['distinct_colors']):
                     return color, color_index
                     
         if color_index >= len(default_colors):
@@ -147,7 +149,6 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
             contour_levels = None
             
         target_styles = super().resolve_target_styles(contour_levels=contour_levels, **kwargs)
-        
         default_colors = self.get_colors()
         color_index = 0
         
@@ -328,6 +329,7 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
     ) -> None:
         """Draw best fit point (minimum likelihood)."""
         styles = mp.concat((self.styles.get('bestfit'), styles), copy=True)
+        styles.setdefault('color', self.get_target_color('bestfit', domain))
         bestfit_x, bestfit_y, bestfit_z = self.get_bestfit(x, y, z)
         
         bestfit_label_fmt = self.get_target_label('bestfit', domain)
@@ -564,13 +566,8 @@ class Likelihood2DPlot(LikelihoodMixin, General2DPlot):
         self.draw_axis_components(ax, xlabel=xlabel, ylabel=ylabel, title=title)
         self.set_axis_range(ax, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
         self.finalize(ax)
-        
-        if legend_order is not None:
-            self.legend_order = legend_order
-        else:
-            self.legend_order = self.get_labelled_legend_domains()
             
         if self.config['draw_legend']:
-            self.draw_legend(ax)
+            self.draw_legend(ax, targets=legend_order)
             
         return ax
\ No newline at end of file
diff --git a/quickstats/plots/likelihood_mixin.py b/quickstats/plots/likelihood_mixin.py
index 6aa8cf27b6c7fa091cc3b782fcf764003762e4a7..bc4872fbbe90b3df8a320a07f953cd321748546f 100644
--- a/quickstats/plots/likelihood_mixin.py
+++ b/quickstats/plots/likelihood_mixin.py
@@ -2,7 +2,7 @@ from __future__ import annotations
 
 from typing import Optional, Dict, Union, Tuple
 
-from quickstats.maths.numerics import str_encode_value
+from quickstats.maths.numerics import str_encode_value, is_integer
 from quickstats.maths.statistics import sigma_to_chi2, confidence_level_to_chi2
 from .abstract_plot import AbstractPlot
 
@@ -42,7 +42,10 @@ class LikelihoodMixin(AbstractPlot):
     def get_level_key(self, level: Union[int, float], use_sigma: bool = False) -> str:
         """Get dictionary key for confidence/sigma level."""
         key = 'sigma_level' if use_sigma else 'confidence_level'
-        level_str = str_encode_value(level)
+        if is_integer(level):
+            level_str = str(int(level))
+        else:
+            level_str = str_encode_value(level)
         return self.config.get('level_key', {}).get(key, '').format(level_str=level_str)
 
     def get_level_specs(
diff --git a/quickstats/plots/template.py b/quickstats/plots/template.py
index 14d2731d2232291b17ade7bc9774536581918ca8..795031c0e8e03c53240ba3dc9d3efc316bbc7590 100644
--- a/quickstats/plots/template.py
+++ b/quickstats/plots/template.py
@@ -17,7 +17,7 @@ import matplotlib.colors as mcolors
 from matplotlib.axes import Axes 
 from matplotlib.axis import Axis
 from matplotlib.artist import Artist
-from matplotlib.patches import Patch, Rectangle, Polygon
+from matplotlib.patches import Patch, Rectangle, Polygon, Ellipse
 from matplotlib.lines import Line2D
 from matplotlib.container import (
     Container,
@@ -1359,7 +1359,6 @@ def remake_handles(
                     **line_styles
                 )
             new_subhandles.append(subhandle)
-            
             if fill_border and isinstance(subhandle, (PolyCollection, BarContainer)):
                 border_style = border_styles or {}
                 border_handle = Rectangle(
@@ -1617,4 +1616,34 @@ def contour_to_shapes(
     for path in contour.get_paths():
         shape = alphashape(path.vertices, alpha)
         shapes.append(shape)
-    return shapes
\ No newline at end of file
+    return shapes
+
+def ellipse_to_shape(
+    ellipse_patch: "Ellipse",
+    num_points: int = 1000
+):
+    from quickstats.core.modules import require_module
+    require_module("shapely")
+    from shapely.geometry import Polygon
+    
+    t = np.linspace(0, 2 * np.pi, num_points)
+    
+    # Ellipse parameters
+    width, height = ellipse_patch.width, ellipse_patch.height
+    angle_rad = np.deg2rad(ellipse_patch.angle)
+    cx, cy = ellipse_patch.center
+    
+    # Parametric equation of ellipse before rotation
+    x = (width / 2) * np.cos(t)
+    y = (height / 2) * np.sin(t)
+    
+    # Rotation matrix
+    R = np.array([
+        [np.cos(angle_rad), -np.sin(angle_rad)],
+        [np.sin(angle_rad),  np.cos(angle_rad)]
+    ])
+    
+    # Apply rotation and translation
+    points = np.column_stack((x, y)) @ R.T + np.array([cx, cy])
+
+    return Polygon(points)
\ No newline at end of file
diff --git a/quickstats/plots/template_styles.py b/quickstats/plots/template_styles.py
index 2c5de4c5a65a93ae6f88b5337768ef6f0ea6d32d..2478121b3fc89a6ac7c69e0d4eaaed8da8552ad5 100644
--- a/quickstats/plots/template_styles.py
+++ b/quickstats/plots/template_styles.py
@@ -138,6 +138,8 @@ REGISTRY['default'] = {
     'vline': {
     },
     'line_collection': {
+    },
+    'ellipse': {
     }
 }