QQ-plot and PP-plot

Both QQ-plot and PP-plot are called the probability plot, but they are different. These plots are intended to compare two distributions, usually at least one of them is empirical. It is to graphically tell how good the two distributions fit.

Assume the two distributions have the cumulative distribution functions \(F(x)\) and \(G(x)\), the PP-plot is to show \(G(x)\) against \(F(x)\) for varying \(x\). Hence the domain and range of the plot is always from 0 to 1, as we are plotting only the range of the cumulative distribution functions.

QQ-plot, however, is to plot the inverse cumulative distribution function \(G^{-1}(x)\) against \(F^{-1}(x)\) for varying \(x\in[0,1]\). Therefore the domain and range of the plot are the support of the cumulative distribution functions \(F(x)\) and \(G(x)\). If we consider the data are empirical, we can see this as the plot of the order statistic of \(G(x)\) against that of \(F(x)\).

Tools

In Python, it is a surprise that matplotlib does not support making PP-plot nor QQ-plot out of the box. However, it should not be difficult to see that we can make use of the order statistics to do the QQ-plot:

import numpy as np
import matplotlib.pyplot as plt

N = 1000 # number of samples
rv_norm = np.random.randn(N) * 2 + 1  # normal with mean 1 s.d. 2
rv_uni = np.random.rand(N) * 8 - 4    # uniform in [-4,4]

plt.scatter(np.sort(rv_uni). np.sort(rv_norm), alpha=0.2, s=2)

PP-plot is a bit more complicated. We need to use interpolation function to achieve that. The idea is that, the CDF of an empirical distribution can be constructed using np.sort(rv_norm) and np.linspace(0,1,N). Then with the other distribution, we can look up the value of the CDF by interpolation using np.interp(x0, x, y), which is to return \(y_0 = f(x_0)\) from the provided curve \(y=f(x)\):

plt.scatter(np.linspace(0,1,N), np.interp(np.sort(rv_uni), np.sort(rv_norm), np.linspace(0,1,N)), alpha=0.2, s=2)

However, there is a fancier tool in Python to do this. statsmodels has functions qqplot() and qqplot_2samples() for doing QQ-plot of one empirical against theoretical normal distribution, and between two empirical distributions, respectively. But it is just a wrapper for the more generic ProbPlot object. For example, this is how we can do the same as above:

import statsmodels.api as sm

_ = sm.ProbPlot(rv_uni) \
      .qqplot(other=sm.ProbPlot(rv_norm), line="r", alpha=0.2, ms=2, lw=1)

plt.show()
_ = sm.ProbPlot(rv_uni) \
      .ppplot(other=sm.ProbPlot(rv_norm), line="r", alpha=0.2, ms=2, lw=1)
plt.show()

its output will go with the regression line (line="r"):

QQ-plot and PP-plot as EDA tool

When we get a table of data the first time, we would like to get some insight from it before further processing it. This is what the exploratory data analysis is about. For a multidimensional data set, my favorite is to run a correlogram to see how the data looks like, visually:

import seaborn as sns
sns.pairplot(df_data)

This graph is generated using Seaborn, a wrapper for matplotlib. We can make the graph prettier, for example, draw the regression line between each pair of data and show the CDF empirically found by KDE instead of histogram:

sns.pairplot(df_data, diag_kind="kde", kind="reg",
             plot_kws={'line_kws':{'color':'red', 'alpha':0.2}, 'scatter_kws': {'alpha': 0.2, 's':4}},
            )
plt.show()

We call the same Seaborn function pairplot() with the kind="reg" (regression) for off-diagonal charts and diag_kind="kde" for on-diagonal charts. This tells you how correlated are any two series and the distribution of each sample. In this graph, we do not see any particularly strong correlation. So how about the series are independent but in the similar distribution? This can be answered by a PP-plot of each pair. Unfortunately, PP-plot and QQ-plot are not supported by Seaborn. Nevertheless we an add this. Here is the patch file, we need only to modify axisgrid.py and regression.py:

diff --git a/seaborn/axisgrid.py b/seaborn/axisgrid.py
index ba70553..7a9d836 100644
--- a/seaborn/axisgrid.py
+++ b/seaborn/axisgrid.py
@@ -1959,6 +1959,7 @@ def pairplot(
     """
     # Avoid circular import
     from .distributions import histplot, kdeplot
+    from .regression import qqplot, ppplot  # Avoid circular import
 
     # Handle deprecations
     if size is not None:
@@ -1992,7 +1993,7 @@ def pairplot(
     # Add the markers here as PairGrid has figured out how many levels of the
     # hue variable are needed and we don't want to duplicate that process
     if markers is not None:
-        if kind == "reg":
+        if kind in ["reg", "pp", "qq"]:
             # Needed until regplot supports style
             if grid.hue_names is None:
                 n_markers = 1
@@ -2020,6 +2021,10 @@ def pairplot(
         diag_kws.setdefault("fill", True)
         diag_kws.setdefault("warn_singular", False)
         grid.map_diag(kdeplot, **diag_kws)
+    elif diag_kind == "pp":
+        grid.map_diag(ppplot, **diag_kws)
+    elif diag_kind == "qq":
+        grid.map_diag(qqplot, **diag_kws)
 
     # Maybe plot on the off-diagonals
     if diag_kind is not None:
@@ -2030,6 +2035,10 @@ def pairplot(
     if kind == "scatter":
         from .relational import scatterplot  # Avoid circular import
         plotter(scatterplot, **plot_kws)
+    elif kind == "qq":
+        plotter(qqplot, **plot_kws)
+    elif kind == "pp":
+        plotter(ppplot, **plot_kws)
     elif kind == "reg":
         from .regression import regplot  # Avoid circular import
         plotter(regplot, **plot_kws)
diff --git a/seaborn/regression.py b/seaborn/regression.py
index ce21927..cc366d1 100644
--- a/seaborn/regression.py
+++ b/seaborn/regression.py
@@ -20,7 +20,7 @@ from .axisgrid import FacetGrid, _facet_docs
 from ._decorators import _deprecate_positional_args
 
 
-__all__ = ["lmplot", "regplot", "residplot"]
+__all__ = ["lmplot", "regplot", "residplot", "ppplot", "qqplot"]
 
 
 class _LinearPlotter(object):
@@ -833,6 +833,91 @@ lmplot.__doc__ = dedent("""\
     """).format(**_regression_docs)
 
 
+@_deprecate_positional_args
+def qqplot(
+    *,
+    x=None, y=None,
+    data=None,
+    x_estimator=None, x_bins=None, x_ci="ci",
+    scatter=True, fit_reg=True, ci=95, n_boot=1000, units=None,
+    seed=None, order=1, logistic=False, lowess=False, robust=False,
+    logx=False, x_partial=None, y_partial=None,
+    truncate=True, dropna=True, x_jitter=None, y_jitter=None,
+    label=None, color=None, marker="o",
+    scatter_kws=None, line_kws=None, ax=None,
+    legend=None
+):
+
+    plotter = _RegressionPlotter(x, y, data, x_estimator, x_bins, x_ci,
+                                 scatter, fit_reg, ci, n_boot, units, seed,
+                                 order, logistic, lowess, robust, logx,
+                                 x_partial, y_partial, truncate, dropna,
+                                 x_jitter, y_jitter, color, label)
+
+    # Manipulate input data for plotting
+    if plotter.x is None:
+        err = "missing x or y in plot data"
+        raise ValueError(err)
+    if plotter.y is None:
+        # set it to normal distribution scaled according to x
+        from scipy.stats import norm
+        plotter.y = norm.ppf(np.linspace(0,1,len(plotter.x)+2)[1:-1])
+        plotter.y = plotter.y * plotter.x.std() + plotter.x.mean()
+    plotter.x = np.sort(plotter.x)
+    plotter.y = np.sort(plotter.y)
+
+    if ax is None:
+        ax = plt.gca()
+
+    scatter_kws = {} if scatter_kws is None else copy.copy(scatter_kws)
+    scatter_kws["marker"] = marker
+    line_kws = {} if line_kws is None else copy.copy(line_kws)
+    plotter.plot(ax, scatter_kws, line_kws)
+    return ax
+
+
+@_deprecate_positional_args
+def ppplot(
+    *,
+    x=None, y=None,
+    data=None,
+    x_estimator=None, x_bins=None, x_ci="ci",
+    scatter=True, fit_reg=True, ci=95, n_boot=1000, units=None,
+    seed=None, order=1, logistic=False, lowess=False, robust=False,
+    logx=False, x_partial=None, y_partial=None,
+    truncate=True, dropna=True, x_jitter=None, y_jitter=None,
+    label=None, color=None, marker="o",
+    scatter_kws=None, line_kws=None, ax=None,
+    legend=None
+):
+
+    plotter = _RegressionPlotter(x, y, data, x_estimator, x_bins, x_ci,
+                                 scatter, fit_reg, ci, n_boot, units, seed,
+                                 order, logistic, lowess, robust, logx,
+                                 x_partial, y_partial, truncate, dropna,
+                                 x_jitter, y_jitter, color, label)
+
+    # Manipulate input data for plotting
+    if plotter.x is None:
+        err = "missing x in plot data"
+        raise ValueError(err)
+    if plotter.y is None:
+        # set it to normal distribution
+        from scipy.stats import norm
+        plotter.y = norm.ppf(np.linspace(0,1,len(plotter.x)+2)[1:-1])
+    linspace = np.linspace(0,1,len(plotter.x))
+    plotter.y = np.interp(np.sort(plotter.x), np.sort(plotter.y), linspace)
+    plotter.x = linspace
+    if plotter.fit_reg:
+        plotter.x_range = (0, 1)
+
+    if ax is None:
+        ax = plt.gca()
+
+    scatter_kws = {} if scatter_kws is None else copy.copy(scatter_kws)
+    scatter_kws["marker"] = marker
+    line_kws = {} if line_kws is None else copy.copy(line_kws)
+    plotter.plot(ax, scatter_kws, line_kws)
+    return ax
+
+
 @_deprecate_positional_args
 def regplot(
     *,

The key changes are the functions sns.ppplot() and sns.qqplot() defined in regression.py, which are modified from the function regplot(). The function regplot() is to do a scatter plot, then make a regression line on top of it. As we saw, PP-plot and QQ-plot are only the modified scatter plot. Therefore we manipulate the data in the plotter using np.sort() and np.interp() before we invoked its plot() function. At this point, these two functions can compare two empirical distributions. However, in the pairplot(), the diagonal charts are handled differently – be it a KDE plot or histogram plot. We can indeed make the PP-plot and QQ-plot a single distribution plot by plotting it against a theoretical normal distribution. The way we do it is to generate one if there is not the second distribution (y): Using the inverse normal CDF function norm.ppf() from scipy, we look for the evenly distributed values from 0 to 1 (clipped the two ends as we know they will be infinite). In case of QQ-plot, the data should be scaled according to the input data to match the mean and standard deviation. The plot will be as follows:

sns.pairplot(df_data, diag_kind="pp", kind="pp",
             plot_kws={'line_kws':{'color':'red', 'alpha':0.2}, 'scatter_kws': {'alpha': 0.2, 's':4}},
             diag_kws={'line_kws':{'color':'red', 'alpha':0.2}, 'scatter_kws': {'alpha': 0.2, 's':4}})
plt.show()

sns.pairplot(df_data, diag_kind="qq", kind="qq",
             plot_kws={'line_kws':{'color':'red', 'alpha':0.2}, 'scatter_kws': {'alpha': 0.2, 's':4}},
             diag_kws={'line_kws':{'color':'red', 'alpha':0.2}, 'scatter_kws': {'alpha': 0.2, 's':4}})
plt.show()

Of course, if you just want one PP-plot, you can use sns.ppplot() directly.

Better solution by PairGrid

The pairplot() function in seaborn is built with the PairGrid object. We can use it directly to make the pair plots but it would not be a single function to achieve that. However, there would be more flexibility. This is an example from the seaborn documentation that allows you to show different chart above and below the diagonal:

# Fancier correlogram
g = sns.PairGrid(df)
g.map_upper(sns.scatterplot)
g.map_lower(sns.kdeplot)
g.map_diag(sns.kdeplot)
plt.show()

There can be PairGrid.map() for all charts in the grid, map_diag() and map_offdiag() for on- and off-diagonal charts, or map_upper() and map_lower() for the upper- and lower-triangular charts of the grid. We can also provide a custom function to the family of PairGrid.map(). The exception of map_diag() is to expect the custom function to take the \(x\) and \(y\) data as positional arguments, while map_diag() takes only \(x\) data. With this, we can define a few functions to be used:

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy.stats import norm

def qqplot(x, y=None, **kwargs):
    """Create Q-Q plot, for use with PairGrid.map()
    Args:
        x, y: Passed as positional arguments, will be pandas series
        kwargs: Can be "color" (tuple of 3 floats) and "label" (str), or other
                keyword arguments passed by user
    """
    # Find scatter locations
    _, xs = stats.probplot(x, fit=False)
    if y is None:
        # set it to normal distribution scaled according to x
        ys = norm.ppf(np.linspace(0,1,len(xs)+2)[1:-1])
        ys = ys * xs.std() + xs.mean()
    else:
        _, ys = stats.probplot(y, fit=False)
    # Find regression line
    p = np.polyfit(xs, ys, 1)
    xr = np.linspace(np.min(xs), np.max(xs), 2)
    yr = np.polyval(p, xr)
    # Plot options
    line_kws = kwargs.pop('line_kws', dict(c="r", ls="--", alpha=0.5))
    scatter_kws = kwargs.pop('scatter_kws', {})
    scatter_kws.update(kwargs) # all remaining keyword arg assumed for scatter
    # Plot
    plt.axline(xy1=(xr[0],yr[0]), xy2=(xr[1],yr[1]), **line_kws)
    sns.scatterplot(x=xs, y=ys, **scatter_kws)

def ppplot(x, y=None, **kwargs):
    """Create P-P plot, for use with PairGrid.map()
    Args:
        x, y: Passed as positional arguments, will be pandas series
        kwargs: Can be "color" (tuple of 3 floats) and "label" (str), or other
                keyword arguments passed by user
    """
    # Find scatter locations
    xs = np.linspace(0, 1, len(x))
    if y is None:
        # set it to normal distribution scaled according to x
        y = norm.ppf(np.linspace(0,1,len(xs)+2)[1:-1])
    ys = np.interp(np.sort(x), np.sort(y), xs)
    # Find regression line
    p = np.polyfit(xs, ys, 1)
    xr = np.linspace(np.min(xs), np.max(xs), 2)
    yr = np.polyval(p, xr)
    # Plot options
    line_kws = kwargs.pop('line_kws', dict(c="r", ls="--", alpha=0.5))
    scatter_kws = kwargs.pop('scatter_kws', {})
    scatter_kws.update(kwargs) # all remaining keyword arg assumed for scatter
    # Plot
    plt.axline(xy1=(xr[0],yr[0]), xy2=(xr[1],yr[1]), **line_kws)
    sns.scatterplot(x=xs, y=ys, **scatter_kws)

def tukeyplot(x, y, **kwargs):
    """Create Tukey mean-difference plot, for use with PairGrid.map()
    Args:
        x, y: Passed as positional arguments, will be pandas series
        kwargs: Can be "color" (tuple of 3 floats) and "label" (str), or other
                keyword arguments passed by user
    """
    # Find locations
    xs = (x+y)/2
    ys = x-y
    mean = np.mean(ys)
    sd = np.std(ys)
    # Plot options
    line_kws = kwargs.pop('line_kws', dict(ls="--", alpha=0.5))
    scatter_kws = kwargs.pop('scatter_kws', {})
    scatter_kws.update(kwargs) # all remaining keyword arg assumed for scatter
    # Plot the mean and 95% ci lines
    xl = np.linspace(np.min(xs), np.max(xs), 4)
    plt.axhline(y=mean, c="b", **line_kws)
    plt.axhline(y=mean+1.96*sd, color="r", **line_kws)
    plt.axhline(y=mean-1.96*sd, color="r", **line_kws)
    sns.scatterplot(x=xs, y=ys, **scatter_kws)

and we can produce a PP-plot with them:

g = sns.PairGrid(df)
g.map_diag(ppplot, line_kws=dict(c="r", alpha=0.3, ls="--"), scatter_kws=dict(alpha=0.5, s=5))
g.map_offdiag(ppplot, scatter_kws=dict(alpha=0.5, s=5))
g.tight_layout()
plt.show()

and a QQ-plot:

g = sns.PairGrid(df)
g.map_diag(qqplot, scatter_kws=dict(alpha=0.5, s=5))
g.map_offdiag(qqplot, scatter_kws=dict(alpha=0.5, s=5))
g.tight_layout()
plt.show()

and also Tukey mean-difference plot with PP-plot on diagonal, but we need to detach the shared axis for the diagonal plots due to different scale:

from itertools import product

g = sns.PairGrid(df)
g.map_offdiag(tukeyplot, scatter_kws=dict(alpha=0.5, s=5))
# Remove shared axis only for diagonal
for i,j,k in product(range(len(df.columns)), repeat=3):
    g.axes[i, i].get_shared_x_axes().remove(g.axes[j, k])
    g.axes[i, i].get_shared_y_axes().remove(g.axes[j, k])
g.map_diag(ppplot, line_kws=dict(alpha=0.3))
g.tight_layout()
plt.show()