diff --git a/README.md b/README.md
index 31d487d5..89fe4090 100644
--- a/README.md
+++ b/README.md
@@ -58,6 +58,7 @@ Hist currently provides everything boost-histogram provides, and the following e
 - Quick plotting routines encourage exploration:
   - `.plot()` provides 1D and 2D plots
   - `.plot2d_full()` shows 1D projects around a 2D plot
+  - `.plot_ratio(...)` make a ratio plot between the histogram and another histogram or callable
   - `.plot_pull(...)` performs a pull plot
   - `.plot_pie()` makes a pie plot
   - `.show()` provides a nice str printout using Histoprint
diff --git a/docs/changelog.rst b/docs/changelog.rst
index c42f2204..a8ba2a3a 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -4,6 +4,10 @@ Changelog
 IN PROGRESS
 --------------------
 
+* Add ``plot_ratio`` to the public API, which allows for making ratio plots between the
+  histogram and either another histogram or a callable.
+  `#161 <https://github.com/scikit-hep/hist/pull/161>`_
+
 * Add frequentist coverage interval support in the ``intervals`` module.
   `#176 <https://github.com/scikit-hep/hist/pull/176>`_
 
diff --git a/src/hist/basehist.py b/src/hist/basehist.py
index 179921f4..cf8d0511 100644
--- a/src/hist/basehist.py
+++ b/src/hist/basehist.py
@@ -35,6 +35,8 @@
     import matplotlib.axes
     from mplhep.plot import Hist1DArtists, Hist2DArtists
 
+    from .plot import FitResultArtists, MainAxisArtists, RatiolikeArtists
+
 InnerIndexing = Union[
     SupportsIndex, str, Callable[[bh.axis.Axis], int], slice, "ellipsis"
 ]
@@ -401,20 +403,45 @@ def plot2d_full(
 
         return hist.plot.plot2d_full(self, ax_dict=ax_dict, **kwargs)
 
+    def plot_ratio(
+        self,
+        other: Union["hist.BaseHist", Callable[[np.ndarray], np.ndarray], str],
+        *,
+        ax_dict: "Optional[Dict[str, matplotlib.axes.Axes]]" = None,
+        **kwargs: Any,
+    ) -> "Tuple[MainAxisArtists, RatiolikeArtists]":
+        """
+        ``plot_ratio`` method for ``BaseHist`` object.
+
+        Return a tuple of artists following a structure of
+        ``(main_ax_artists, subplot_ax_artists)``
+        """
+
+        import hist.plot
+
+        return hist.plot._plot_ratiolike(
+            self, other, ax_dict=ax_dict, view="ratio", **kwargs
+        )
+
     def plot_pull(
         self,
-        func: Callable[[np.ndarray], np.ndarray],
+        func: Union[Callable[[np.ndarray], np.ndarray], str],
         *,
         ax_dict: "Optional[Dict[str, matplotlib.axes.Axes]]" = None,
         **kwargs: Any,
-    ) -> "Tuple[matplotlib.axes.Axes, matplotlib.axes.Axes]":
+    ) -> "Tuple[FitResultArtists, RatiolikeArtists]":
         """
-        Plot_pull method for BaseHist object.
+        ``plot_pull`` method for ``BaseHist`` object.
+
+        Return a tuple of artists following a structure of
+        ``(main_ax_artists, subplot_ax_artists)``
         """
 
         import hist.plot
 
-        return hist.plot.plot_pull(self, func, ax_dict=ax_dict, **kwargs)
+        return hist.plot._plot_ratiolike(
+            self, func, ax_dict=ax_dict, view="pull", **kwargs
+        )
 
     def plot_pie(
         self,
diff --git a/src/hist/plot.py b/src/hist/plot.py
index 2cec5e03..d69d1355 100644
--- a/src/hist/plot.py
+++ b/src/hist/plot.py
@@ -1,11 +1,25 @@
 import inspect
 import sys
-from typing import Any, Callable, Dict, Iterable, Optional, Set, Tuple, Union
+from typing import (
+    Any,
+    Callable,
+    Dict,
+    Iterable,
+    List,
+    NamedTuple,
+    Optional,
+    Set,
+    Tuple,
+    Union,
+)
 
 import numpy as np
 
 import hist
 
+from .intervals import ratio_uncertainty
+from .typing import Literal
+
 try:
     import matplotlib.axes
     import matplotlib.patches as patches
@@ -19,8 +33,42 @@
     )
     raise
 
+__all__ = (
+    "histplot",
+    "hist2dplot",
+    "plot2d_full",
+    "plot_ratio_array",
+    "plot_pull_array",
+    "plot_pie",
+)
+
+
+class FitResultArtists(NamedTuple):
+    line: matplotlib.lines.Line2D
+    errorbar: matplotlib.container.ErrorbarContainer
+    band: matplotlib.collections.PolyCollection
+
+
+class RatioErrorbarArtists(NamedTuple):
+    line: matplotlib.lines.Line2D
+    errorbar: matplotlib.container.ErrorbarContainer
+
+
+class RatioBarArtists(NamedTuple):
+    line: matplotlib.lines.Line2D
+    dots: matplotlib.collections.PathCollection
+    bar: matplotlib.container.BarContainer
+
+
+class PullArtists(NamedTuple):
+    bar: matplotlib.container.BarContainer
+    patch_artist: List[matplotlib.patches.Rectangle]
+
 
-__all__ = ("histplot", "hist2dplot", "plot2d_full", "plot_pull", "plot_pie")
+MainAxisArtists = Union[FitResultArtists, Hist1DArtists]
+
+RatioArtists = Union[RatioErrorbarArtists, RatioBarArtists]
+RatiolikeArtists = Union[RatioArtists, PullArtists]
 
 
 def __dir__() -> Tuple[str, ...]:
@@ -134,9 +182,9 @@ def fnll(v: Iterable[np.ndarray]) -> float:
 def plot2d_full(
     self: hist.BaseHist,
     *,
-    ax_dict: "Optional[Dict[str, matplotlib.axes.Axes]]" = None,
+    ax_dict: Optional[Dict[str, matplotlib.axes.Axes]] = None,
     **kwargs: Any,
-) -> "Tuple[Hist2DArtists, Hist1DArtists, Hist1DArtists]":
+) -> Tuple[Hist2DArtists, Hist1DArtists, Hist1DArtists]:
     """
     Plot2d_full method for BaseHist object.
 
@@ -214,19 +262,261 @@ def plot2d_full(
     return main_art, top_art, side_art
 
 
-def plot_pull(
+def _construct_gaussian_callable(
+    __hist: hist.BaseHist,
+) -> Callable[[np.ndarray], np.ndarray]:
+    x_values = __hist.axes[0].centers
+    hist_values = __hist.values()
+
+    # gaussian with reasonable initial guesses for parameters
+    constant = float(hist_values.max())
+    mean = (hist_values * x_values).sum() / hist_values.sum()
+    sigma = (hist_values * np.square(x_values - mean)).sum() / hist_values.sum()
+
+    # gauss is a closure that will get evaluated in _fit_callable_to_hist
+    def gauss(
+        x: np.ndarray,
+        constant: float = constant,
+        mean: float = mean,
+        sigma: float = sigma,
+    ) -> np.ndarray:
+        # Note: Force np.ndarray type as numpy ufuncs have type "Any"
+        ret: np.ndarray = constant * np.exp(
+            -np.square(x - mean) / (2 * np.square(sigma))
+        )
+        return ret
+
+    return gauss
+
+
+def _fit_callable_to_hist(
+    model: Callable[[np.ndarray], np.ndarray],
+    histogram: hist.BaseHist,
+    likelihood: bool = False,
+) -> "Tuple[np.ndarray, np.ndarray, np.ndarray, Tuple[Tuple[float, ...], np.ndarray]]":
+    """
+    Fit a model, a callable function, to the histogram values.
+    """
+    variances = histogram.variances()
+    if variances is None:
+        raise RuntimeError(
+            "Cannot compute from a variance-less histogram, try a Weight storage"
+        )
+    hist_uncert = np.sqrt(variances)
+
+    # Infer best fit model parameters and covariance matrix
+    xdata = histogram.axes[0].centers
+    popt, pcov = _curve_fit_wrapper(
+        model, xdata, histogram.values(), hist_uncert, likelihood=likelihood
+    )
+    model_values = model(xdata, *popt)
+
+    if np.isfinite(pcov).all():
+        n_samples = 100
+        vopts = np.random.multivariate_normal(popt, pcov, n_samples)
+        sampled_ydata = np.vstack([model(xdata, *vopt).T for vopt in vopts])
+        model_uncert = np.nanstd(sampled_ydata, axis=0)
+    else:
+        model_uncert = np.zeros_like(hist_uncert)
+
+    return model_values, model_uncert, hist_uncert, (popt, pcov)
+
+
+def _plot_fit_result(
+    __hist: hist.BaseHist,
+    model_values: np.ndarray,
+    model_uncert: np.ndarray,
+    ax: matplotlib.axes.Axes,
+    eb_kwargs: Dict[str, Any],
+    fp_kwargs: Dict[str, Any],
+    ub_kwargs: Dict[str, Any],
+) -> FitResultArtists:
+    """
+    Plot fit of model to histogram data
+    """
+    x_values = __hist.axes[0].centers
+    variances = __hist.variances()
+    if variances is None:
+        raise RuntimeError(
+            "Cannot compute from a variance-less histogram, try a Weight storage"
+        )
+    hist_uncert = np.sqrt(variances)
+
+    errorbars = ax.errorbar(x_values, __hist.values(), hist_uncert, **eb_kwargs)
+
+    # Ensure zorder draws data points above model
+    line_zorder = errorbars[0].get_zorder() - 1
+    (line,) = ax.plot(x_values, model_values, **fp_kwargs, zorder=line_zorder)
+
+    # Uncertainty band for fitted function
+    ub_kwargs.setdefault("color", line.get_color())
+    if ub_kwargs["color"] == line.get_color():
+        ub_kwargs.setdefault("alpha", 0.3)
+    uncertainty_band = ax.fill_between(
+        x_values,
+        model_values - model_uncert,
+        model_values + model_uncert,
+        **ub_kwargs,
+    )
+
+    return FitResultArtists(line, errorbars, uncertainty_band)
+
+
+def plot_ratio_array(
+    __hist: hist.BaseHist,
+    ratio: np.ndarray,
+    ratio_uncert: np.ndarray,
+    ax: matplotlib.axes.Axes,
+    **kwargs: Any,
+) -> RatioArtists:
+    """
+    Plot a ratio plot on the given axes
+    """
+    x_values = __hist.axes[0].centers
+    left_edge = __hist.axes.edges[0][0]
+    right_edge = __hist.axes.edges[-1][-1]
+
+    # Set 0 and inf to nan to hide during plotting
+    ratio[ratio == 0] = np.nan
+    ratio[np.isinf(ratio)] = np.nan
+
+    central_value = kwargs.pop("central_value", 1.0)
+    central_value_artist = ax.axhline(
+        central_value, color="black", linestyle="dashed", linewidth=1.0
+    )
+
+    # Type now due to control flow
+    axis_artists: Union[RatioErrorbarArtists, RatioBarArtists]
+
+    uncert_draw_type = kwargs.pop("uncert_draw_type", "line")
+    if uncert_draw_type == "line":
+        errorbar_artists = ax.errorbar(
+            x_values,
+            ratio,
+            yerr=ratio_uncert,
+            color="black",
+            marker="o",
+            linestyle="none",
+        )
+        axis_artists = RatioErrorbarArtists(central_value_artist, errorbar_artists)
+    elif uncert_draw_type == "bar":
+        bar_width = (right_edge - left_edge) / len(ratio)
+
+        bar_top = ratio + ratio_uncert[1]
+        bar_bottom = ratio - ratio_uncert[0]
+        # bottom can't be nan
+        bar_bottom[np.isnan(bar_bottom)] = 0
+        bar_height = bar_top - bar_bottom
+
+        _ratio_points = ax.scatter(x_values, ratio, color="black")
+
+        # Ensure zorder draws data points above uncertainty bars
+        bar_zorder = _ratio_points.get_zorder() - 1
+        bar_artists = ax.bar(
+            x_values,
+            height=bar_height,
+            width=bar_width,
+            bottom=bar_bottom,
+            fill=False,
+            linewidth=0,
+            edgecolor="gray",
+            hatch=3 * "/",
+            zorder=bar_zorder,
+        )
+        axis_artists = RatioBarArtists(central_value_artist, _ratio_points, bar_artists)
+
+    ratio_ylim = kwargs.pop("ylim", None)
+    if ratio_ylim is None:
+        # plot centered around central value with a scaled view range
+        # the value _with_ the uncertainty in view is important so base
+        # view range on extrema of value +/- uncertainty
+        valid_ratios_idx = np.where(np.isnan(ratio) == False)  # noqa: E712
+        valid_ratios = ratio[valid_ratios_idx]
+        extrema = np.array(
+            [
+                valid_ratios - ratio_uncert[0][valid_ratios_idx],
+                valid_ratios + ratio_uncert[1][valid_ratios_idx],
+            ]
+        )
+        max_delta = np.max(np.abs(extrema - central_value))
+        ratio_extrema = np.abs(max_delta + central_value)
+
+        _alpha = 2.0
+        scaled_offset = max_delta + (max_delta / (_alpha * ratio_extrema))
+        ratio_ylim = [central_value - scaled_offset, central_value + scaled_offset]
+
+    ax.set_xlim(left_edge, right_edge)
+    ax.set_ylim(bottom=ratio_ylim[0], top=ratio_ylim[1])
+
+    ax.set_xlabel(__hist.axes[0].label)
+    ax.set_ylabel(kwargs.pop("ylabel", "Ratio"))
+
+    return axis_artists
+
+
+def plot_pull_array(
+    __hist: hist.BaseHist,
+    pulls: np.ndarray,
+    ax: matplotlib.axes.Axes,
+    bar_kwargs: Dict[str, Any],
+    pp_kwargs: Dict[str, Any],
+) -> PullArtists:
+    """
+    Plot a pull plot on the given axes
+    """
+    x_values = __hist.axes[0].centers
+    left_edge = __hist.axes.edges[0][0]
+    right_edge = __hist.axes.edges[-1][-1]
+
+    # Pull: plot the pulls using Matplotlib bar method
+    width = (right_edge - left_edge) / len(pulls)
+    bar_artists = ax.bar(x_values, pulls, width=width, **bar_kwargs)
+
+    pp_num = pp_kwargs.pop("num", 5)
+    patch_height = max(np.abs(pulls)) / pp_num
+    patch_width = width * len(pulls)
+    patch_artists = []
+    for i in range(pp_num):
+        # gradient color patches
+        if "alpha" in pp_kwargs:
+            pp_kwargs["alpha"] *= np.power(0.618, i)
+        else:
+            pp_kwargs["alpha"] = 0.5 * np.power(0.618, i)
+
+        upRect_startpoint = (left_edge, i * patch_height)
+        upRect = patches.Rectangle(
+            upRect_startpoint, patch_width, patch_height, **pp_kwargs
+        )
+        ax.add_patch(upRect)
+        downRect_startpoint = (left_edge, -(i + 1) * patch_height)
+        downRect = patches.Rectangle(
+            downRect_startpoint, patch_width, patch_height, **pp_kwargs
+        )
+        ax.add_patch(downRect)
+        patch_artists.append((downRect, upRect))
+
+    ax.set_xlim(left_edge, right_edge)
+
+    ax.set_xlabel(__hist.axes[0].label)
+    ax.set_ylabel("Pull")
+
+    return PullArtists(bar_artists, patch_artists)
+
+
+def _plot_ratiolike(
     self: hist.BaseHist,
-    func: Union[Callable[[np.ndarray], np.ndarray], str],
+    other: Union[hist.BaseHist, Callable[[np.ndarray], np.ndarray], str],
     likelihood: bool = False,
     *,
-    ax_dict: "Optional[Dict[str, matplotlib.axes.Axes]]" = None,
+    ax_dict: Optional[Dict[str, matplotlib.axes.Axes]] = None,
+    view: Literal["ratio", "pull"],
     fit_fmt: Optional[str] = None,
     **kwargs: Any,
-) -> "Tuple[matplotlib.axes.Axes, matplotlib.axes.Axes]":
+) -> Tuple[MainAxisArtists, RatiolikeArtists]:
     r"""
-    Plot_pull method for BaseHist object.
+    Plot ratio-like plots (ratio plots and pull plots) for BaseHist
 
-    fit_fmt can be a string such as r"{name} = {value:.3g} $\pm$ {error:.3g}"
+    ``fit_fmt`` can be a string such as ``r"{name} = {value:.3g} $\pm$ {error:.3g}"``
     """
 
     try:
@@ -234,23 +524,24 @@ def plot_pull(
         from scipy.optimize import curve_fit  # noqa: F401
     except ModuleNotFoundError:
         print(
-            "Hist.plot_pull requires scipy and iminuit. Please install hist[plot] or manually install dependencies.",
+            f"Hist.plot_{view} requires scipy and iminuit. Please install hist[plot] or manually install dependencies.",
             file=sys.stderr,
         )
         raise
 
-    # Type judgement
-    if not callable(func) and not type(func) in [str]:
-        msg = f"Parameter func must be callable or a string for {self.__class__.__name__} in plot pull"
-        raise TypeError(msg)
-
     if self.ndim != 1:
-        raise TypeError("Only 1D-histogram supports pull plot, try projecting to 1D")
+        raise TypeError(
+            f"Only 1D-histogram supports ratio plot, try projecting {self.__class__.__name__} to 1D"
+        )
+    if isinstance(other, hist.hist.Hist) and other.ndim != 1:
+        raise TypeError(
+            f"Only 1D-histogram supports ratio plot, try projecting other={other.__class__.__name__} to 1D"
+        )
 
     if ax_dict:
         try:
             main_ax = ax_dict["main_ax"]
-            pull_ax = ax_dict["pull_ax"]
+            subplot_ax = ax_dict[f"{view}_ax"]
         except KeyError:
             raise ValueError("All axes should be all given or none at all")
     else:
@@ -258,141 +549,121 @@ def plot_pull(
         grid = fig.add_gridspec(2, 1, hspace=0, height_ratios=[3, 1])
 
         main_ax = fig.add_subplot(grid[0])
-        pull_ax = fig.add_subplot(grid[1], sharex=main_ax)
+        subplot_ax = fig.add_subplot(grid[1], sharex=main_ax)
         plt.setp(main_ax.get_xticklabels(), visible=False)
 
-    # Computation and Fit
-    xdata = self.axes[0].centers
-    ydata = self.values()
-    variances = self.variances()
-    if variances is None:
-        raise RuntimeError(
-            "Cannot compute from a variance-less histogram, try a Weight storage"
-        )
-    yerr = np.sqrt(variances)
-
-    if isinstance(func, str):
-        if func in {"gauss", "gaus"}:
-            # gaussian with reasonable initial guesses for parameters
-            constant = float(ydata.max())
-            mean = (ydata * xdata).sum() / ydata.sum()
-            sigma = (ydata * (xdata - mean) ** 2.0).sum() / ydata.sum()
-
-            def func(
-                x: np.ndarray,
-                constant: float = constant,
-                mean: float = mean,
-                sigma: float = sigma,
-            ) -> np.ndarray:
-                return constant * np.exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))  # type: ignore
-
-        else:
-            func = _expr_to_lambda(func)
-
-    assert not isinstance(func, str)
-
-    parnames = list(inspect.signature(func).parameters)[1:]
-
-    # Compute fit values: using func as fit model
-    popt, pcov = _curve_fit_wrapper(func, xdata, ydata, yerr, likelihood=likelihood)
-    perr = np.diag(pcov) ** 0.5
-    yfit = func(self.axes[0].centers, *popt)
-
-    if np.isfinite(pcov).all():
-        nsamples = 100
-        vopts = np.random.multivariate_normal(popt, pcov, nsamples)
-        sampled_ydata = np.vstack([func(xdata, *vopt).T for vopt in vopts])
-        yfiterr = np.nanstd(sampled_ydata, axis=0)
-    else:
-        yfiterr = np.zeros_like(yerr)
-
-    # Compute pulls: containing no INF values
-    with np.errstate(divide="ignore"):
-        pulls = (ydata - yfit) / yerr
-
-    pulls[np.isnan(pulls)] = 0
-    pulls[np.isinf(pulls)] = 0
-
     # Keyword Argument Conversion: convert the kwargs to several independent args
-
     # error bar keyword arguments
     eb_kwargs = _filter_dict(kwargs, "eb_")
     eb_kwargs.setdefault("label", "Histogram Data")
+    # Use "fmt" over "marker" to avoid UserWarning on keyword precedence
     eb_kwargs.setdefault("fmt", "o")
+    eb_kwargs.setdefault("linestyle", "none")
 
     # fit plot keyword arguments
-    label = "Fit"
-    if fit_fmt is not None:
-        for name, value, error in zip(parnames, popt, perr):
-            label += "\n  "
-            label += fit_fmt.format(name=name, value=value, error=error)
     fp_kwargs = _filter_dict(kwargs, "fp_")
-    fp_kwargs.setdefault("label", label)
+    fp_kwargs.setdefault("label", "Counts")
+
+    # bar plot keyword arguments
+    bar_kwargs = _filter_dict(kwargs, "bar_", ignore={"bar_width"})
 
     # uncertainty band keyword arguments
     ub_kwargs = _filter_dict(kwargs, "ub_")
     ub_kwargs.setdefault("label", "Uncertainty")
-    ub_kwargs.setdefault("alpha", 0.5)
 
-    # bar plot keyword arguments
-    bar_kwargs = _filter_dict(kwargs, "bar_", ignore={"bar_width"})
+    # ratio plot keyword arguments
+    rp_kwargs = _filter_dict(kwargs, "rp_")
+    rp_kwargs.setdefault("uncertainty_type", "poisson")
+    rp_kwargs.setdefault("legend_loc", "best")
+    rp_kwargs.setdefault("num_label", None)
+    rp_kwargs.setdefault("denom_label", None)
 
     # patch plot keyword arguments
-    pp_kwargs = _filter_dict(kwargs, "pp_", ignore={"pp_num"})
-    pp_num = kwargs.pop("pp_num", 5)
+    pp_kwargs = _filter_dict(kwargs, "pp_")
 
     # Judge whether some arguments are left
     if kwargs:
         raise ValueError(f"{set(kwargs)}' not needed")
 
-    # Main: plot the pulls using Matplotlib errorbar and plot methods
-    main_ax.errorbar(self.axes.centers[0], ydata, yerr, **eb_kwargs)
-
-    (line,) = main_ax.plot(self.axes.centers[0], yfit, **fp_kwargs)
-
-    # Uncertainty band
-    ub_kwargs.setdefault("color", line.get_color())
-    main_ax.fill_between(
-        self.axes.centers[0],
-        yfit - yfiterr,
-        yfit + yfiterr,
-        **ub_kwargs,
-    )
-    main_ax.legend(loc=0)
-    main_ax.set_ylabel("Counts")
-
-    # Pull: plot the pulls using Matplotlib bar method
-    left_edge = self.axes.edges[0][0]
-    right_edge = self.axes.edges[-1][-1]
-    width = (right_edge - left_edge) / len(pulls)
-    pull_ax.bar(self.axes.centers[0], pulls, width=width, **bar_kwargs)
+    main_ax.set_ylabel(fp_kwargs["label"])
 
-    patch_height = max(np.abs(pulls)) / pp_num
-    patch_width = width * len(pulls)
-    for i in range(pp_num):
-        # gradient color patches
-        if "alpha" in pp_kwargs:
-            pp_kwargs["alpha"] *= np.power(0.618, i)
+    # Computation and Fit
+    hist_values = self.values()
+
+    main_ax_artists: MainAxisArtists  # Type now due to control flow
+    if callable(other) or isinstance(other, str):
+        if isinstance(other, str):
+            if other in {"gauss", "gaus", "normal"}:
+                other = _construct_gaussian_callable(self)
+            else:
+                other = _expr_to_lambda(other)
+
+        (
+            compare_values,
+            model_uncert,
+            hist_values_uncert,
+            bestfit_result,
+        ) = _fit_callable_to_hist(other, self, likelihood)
+
+        if fit_fmt is not None:
+            parnames = list(inspect.signature(other).parameters)[1:]
+            popt, pcov = bestfit_result
+            perr = np.sqrt(np.diag(pcov))
+
+            fp_label = "Fit"
+            for name, value, error in zip(parnames, popt, perr):
+                fp_label += "\n  "
+                fp_label += fit_fmt.format(name=name, value=value, error=error)
+            fp_kwargs["label"] = fp_label
         else:
-            pp_kwargs["alpha"] = 0.5 * np.power(0.618, i)
-
-        upRect_startpoint = (left_edge, i * patch_height)
-        upRect = patches.Rectangle(
-            upRect_startpoint, patch_width, patch_height, **pp_kwargs
-        )
-        pull_ax.add_patch(upRect)
-        downRect_startpoint = (left_edge, -(i + 1) * patch_height)
-        downRect = patches.Rectangle(
-            downRect_startpoint, patch_width, patch_height, **pp_kwargs
+            fp_kwargs["label"] = "Fitted value"
+
+        main_ax_artists = _plot_fit_result(
+            self,
+            model_values=compare_values,
+            model_uncert=model_uncert,
+            ax=main_ax,
+            eb_kwargs=eb_kwargs,
+            fp_kwargs=fp_kwargs,
+            ub_kwargs=ub_kwargs,
         )
-        pull_ax.add_patch(downRect)
+    else:
+        compare_values = other.values()
+
+        self_artists = histplot(self, ax=main_ax, label=rp_kwargs["num_label"])
+        other_artists = histplot(other, ax=main_ax, label=rp_kwargs["denom_label"])
+
+        main_ax_artists = self_artists, other_artists
+
+    subplot_ax_artists: RatiolikeArtists  # Type now due to control flow
+    # Compute ratios: containing no INF values
+    with np.errstate(divide="ignore", invalid="ignore"):
+        if view == "ratio":
+            ratios = hist_values / compare_values
+            ratio_uncert = ratio_uncertainty(
+                num=hist_values,
+                denom=compare_values,
+                uncertainty_type=rp_kwargs["uncertainty_type"],
+            )
+            # ratio: plot the ratios using Matplotlib errorbar or bar
+            subplot_ax_artists = plot_ratio_array(
+                self, ratios, ratio_uncert, ax=subplot_ax, **rp_kwargs
+            )
+
+        elif view == "pull":
+            pulls = (hist_values - compare_values) / hist_values_uncert
 
-    plt.xlim(left_edge, right_edge)
+            pulls[np.isnan(pulls) | np.isinf(pulls)] = 0
+
+            # Pass dicts instead of unpacking to avoid conflicts
+            subplot_ax_artists = plot_pull_array(
+                self, pulls, ax=subplot_ax, bar_kwargs=bar_kwargs, pp_kwargs=pp_kwargs
+            )
 
-    pull_ax.set_xlabel(self.axes[0].label)
-    pull_ax.set_ylabel("Pull")
+    if main_ax.get_legend_handles_labels()[0]:  # Don't plot an empty legend
+        main_ax.legend(loc=rp_kwargs["legend_loc"])
 
-    return main_ax, pull_ax
+    return main_ax_artists, subplot_ax_artists
 
 
 def get_center(x: Union[str, int, Tuple[float, float]]) -> Union[str, float]:
@@ -405,7 +676,7 @@ def get_center(x: Union[str, int, Tuple[float, float]]) -> Union[str, float]:
 def plot_pie(
     self: hist.BaseHist,
     *,
-    ax: "Optional[matplotlib.axes.Axes]" = None,
+    ax: Optional[matplotlib.axes.Axes] = None,
     **kwargs: Any,
 ) -> Any:
 
diff --git a/tests/baseline/test_image_plot_pull.png b/tests/baseline/test_image_plot_pull.png
index 7bfde203..b5354860 100644
Binary files a/tests/baseline/test_image_plot_pull.png and b/tests/baseline/test_image_plot_pull.png differ
diff --git a/tests/baseline/test_image_plot_ratio_callable.png b/tests/baseline/test_image_plot_ratio_callable.png
new file mode 100644
index 00000000..ac6e15e9
Binary files /dev/null and b/tests/baseline/test_image_plot_ratio_callable.png differ
diff --git a/tests/baseline/test_image_plot_ratio_hist.png b/tests/baseline/test_image_plot_ratio_hist.png
new file mode 100644
index 00000000..ddc83f22
Binary files /dev/null and b/tests/baseline/test_image_plot_ratio_hist.png differ
diff --git a/tests/test_plot.py b/tests/test_plot.py
index 1b6d55c2..5a4d3563 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -534,10 +534,6 @@ def pdf(x, a=1 / np.sqrt(2 * np.pi), x0=0, sigma=1, offset=0):
 
     assert h.plot_pull(pdf_str)
 
-    assert h.plot_pull("gauss")
-
-    assert h.plot_pull("gauss", likelihood=True)
-
     # dimension error
     hh = NamedHist(
         axis.Regular(
@@ -600,6 +596,25 @@ def pdf(x, a=1 / np.sqrt(2 * np.pi), x0=0, sigma=1, offset=0):
     plt.close("all")
 
 
+@pytest.mark.parametrize("str_alias", ["normal", "gauss", "gaus"])
+@pytest.mark.parametrize("use_likelihood", [True, False])
+def test_ratiolike_str_alias(str_alias, use_likelihood):
+    """
+    Test str alias for callable in plot_ratio and plot_pull
+    """
+
+    np.random.seed(42)
+
+    h = NamedHist(
+        axis.Regular(
+            50, -4, 4, name="S", label="s [units]", underflow=False, overflow=False
+        )
+    ).fill(S=np.random.normal(size=10))
+
+    assert h.plot_ratio(str_alias, likelihood=use_likelihood)
+    assert h.plot_pull(str_alias, likelihood=use_likelihood)
+
+
 @pytest.mark.mpl_image_compare(baseline_dir="baseline", savefig_kwargs={"dpi": 50})
 def test_image_plot_pull():
     """
@@ -620,7 +635,69 @@ def pdf(x, a=1 / np.sqrt(2 * np.pi), x0=0, sigma=1, offset=0):
 
     fig = plt.figure()
 
-    assert h.plot_pull(pdf, fit_fmt=r"{name} = {value:.3g} $\pm$ {error:.3g}")
+    assert h.plot_pull(
+        pdf,
+        eb_color="black",
+        fp_color="blue",
+        ub_color="lightblue",
+        fit_fmt=r"{name} = {value:.3g} $\pm$ {error:.3g}",
+    )
+
+    return fig
+
+
+@pytest.mark.mpl_image_compare(baseline_dir="baseline", savefig_kwargs={"dpi": 50})
+def test_image_plot_ratio_hist():
+    """
+    Test plot_pull by comparing against a reference image generated via
+    `pytest --mpl-generate-path=tests/baseline`
+    """
+
+    np.random.seed(42)
+
+    hist_1 = Hist(
+        axis.Regular(
+            50, -5, 5, name="X", label="x [units]", underflow=False, overflow=False
+        )
+    ).fill(np.random.normal(size=1000))
+    hist_2 = Hist(
+        axis.Regular(
+            50, -5, 5, name="X", label="x [units]", underflow=False, overflow=False
+        )
+    ).fill(np.random.normal(size=1700))
+
+    fig = plt.figure()
+
+    assert hist_1.plot_ratio(
+        hist_2, rp_num_label="numerator", rp_denom_label="denominator"
+    )
+
+    return fig
+
+
+@pytest.mark.mpl_image_compare(baseline_dir="baseline", savefig_kwargs={"dpi": 50})
+def test_image_plot_ratio_callable():
+    """
+    Test plot_pull by comparing against a reference image generated via
+    `pytest --mpl-generate-path=tests/baseline`
+    """
+
+    np.random.seed(42)
+
+    hist_1 = Hist(
+        axis.Regular(
+            50, -5, 5, name="X", label="x [units]", underflow=False, overflow=False
+        )
+    ).fill(np.random.normal(size=1000))
+
+    def model(x, a=1 / np.sqrt(2 * np.pi), x0=0, sigma=1, offset=0):
+        return a * np.exp(-((x - x0) ** 2) / (2 * sigma ** 2)) + offset
+
+    fig = plt.figure()
+
+    assert hist_1.plot_ratio(
+        model, eb_color="black", fp_color="blue", ub_color="lightblue"
+    )
 
     return fig