1
$\begingroup$

Suppose I have two samples from unknown distributions, and I want to produce a quantile-quantile plot of one against the other. What's the easiest way to do this in Python?

There is scipy.stats.probplot, and something similar in statsmodels, but both produce a q-q plot of data versus some predefined distribution, which is not what I want.

$\endgroup$
6
$\begingroup$

A simple one-off solution:

import numpy as np
import matplotlib.pyplot as plt

# Setup
rng = np.random.RandomState(0)  # Seed RNG for replicability
n = 100  # Number of samples to draw

# Generate data
x = rng.normal(size=n)  # Sample 1: X ~ N(0, 1)
y = rng.standard_t(df=5, size=n)  # Sample 2: Y ~ t(5)

# Quantile-quantile plot
plt.figure()
plt.scatter(np.sort(x), np.sort(y))
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
plt.close()

which produces

qqplot

Better yet, you can modularize this code by writing a q-q plot function, something like the following:

import numbers

import numpy as np
import matplotlib.pyplot as plt


def qqplot(x, y, quantiles=None, interpolation='nearest', ax=None, rug=False,
           rug_length=0.05, rug_kwargs=None, **kwargs):
    """Draw a quantile-quantile plot for `x` versus `y`.

    Parameters
    ----------
    x, y : array-like
        One-dimensional numeric arrays.

    ax : matplotlib.axes.Axes, optional
        Axes on which to plot. If not provided, the current axes will be used.

    quantiles : int or array-like, optional
        Quantiles to include in the plot. This can be an array of quantiles, in
        which case only the specified quantiles of `x` and `y` will be plotted.
        If this is an int `n`, then the quantiles will be `n` evenly spaced
        points between 0 and 1. If this is None, then `min(len(x), len(y))`
        evenly spaced quantiles between 0 and 1 will be computed.

    interpolation : {‘linear’, ‘lower’, ‘higher’, ‘midpoint’, ‘nearest’}
        Specify the interpolation method used to find quantiles when `quantiles`
        is an int or None. See the documentation for numpy.quantile().

    rug : bool, optional
        If True, draw a rug plot representing both samples on the horizontal and
        vertical axes. If False, no rug plot is drawn.

    rug_length : float in [0, 1], optional
        Specifies the length of the rug plot lines as a fraction of the total
        vertical or horizontal length.

    rug_kwargs : dict of keyword arguments
        Keyword arguments to pass to matplotlib.axes.Axes.axvline() and
        matplotlib.axes.Axes.axhline() when drawing rug plots.

    kwargs : dict of keyword arguments
        Keyword arguments to pass to matplotlib.axes.Axes.scatter() when drawing
        the q-q plot.
    """
    # Get current axes if none are provided
    if ax is None:
        ax = plt.gca()

    if quantiles is None:
        quantiles = min(len(x), len(y))

    # Compute quantiles of the two samples
    if isinstance(quantiles, numbers.Integral):
        quantiles = np.linspace(start=0, stop=1, num=int(quantiles))
    else:
        quantiles = np.atleast_1d(np.sort(quantiles))
    x_quantiles = np.quantile(x, quantiles, interpolation=interpolation)
    y_quantiles = np.quantile(y, quantiles, interpolation=interpolation)

    # Draw the rug plots if requested
    if rug:
        # Default rug plot settings
        rug_x_params = dict(ymin=0, ymax=rug_length, c='gray', alpha=0.5)
        rug_y_params = dict(xmin=0, xmax=rug_length, c='gray', alpha=0.5)

        # Override default setting by any user-specified settings
        if rug_kwargs is not None:
            rug_x_params.update(rug_kwargs)
            rug_y_params.update(rug_kwargs)

        # Draw the rug plots
        for point in x:
            ax.axvline(point, **rug_x_params)
        for point in y:
            ax.axhline(point, **rug_y_params)

    # Draw the q-q plot
    ax.scatter(x_quantiles, y_quantiles, **kwargs)

This can be used as follows:

# Setup
rng = np.random.RandomState(0)  # Seed RNG for replicability

# Example 1: samples of the same length
n = 100  # Number of samples to draw
x = rng.normal(size=n)  # Sample 1: X ~ N(0, 1)
y = rng.standard_t(df=5, size=n)  # Sample 2: Y ~ t(5)

# Draw quantile-quantile plot
plt.figure()
qqplot(x, y, c='r', alpha=0.5, edgecolor='k')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Two samples of the same length')
plt.show()
plt.close()

# Example 2: samples of different lengths
n_u = 50  # Number of U samples to draw
n_v = 100  # Number of V samples to draw
u = rng.normal(size=n_u)  # Sample 1: U ~ N(0, 1)
v = rng.standard_t(df=5, size=n_v)  # Sample 2: V ~ t(5)

# Draw quantile-quantile plot
plt.figure()
qqplot(u, v, c='r', alpha=0.5, edgecolor='k')
plt.xlabel('U')
plt.ylabel('V')
plt.title('Two samples of different lengths')
plt.show()
plt.close()

# Draw quantile-quantile plot with rug plot
plt.figure()
qqplot(u, v, c='r', alpha=0.5, edgecolor='k', rug=True)
plt.xlabel('U')
plt.ylabel('V')
plt.title('Two samples of different lengths, with rug plot')
plt.show()
plt.close()

which produces

$\endgroup$
  • 1
    $\begingroup$ Here is a somewhat more general thing (works for different-length inputs) (feel free to incorporate it into your answer, I don't have enough rep here to edit). def qqplot(x, y, numquant=None, ax=None, **kwargs): if ax is None: ax = plt.gca() if numquant is None: numquant = min(len(x), len(y)) x1=np.quantile(x, np.linspace(0, 1, numquant)) y1=np.quantile(y, np.linspace(0, 1, numquant)) ax.scatter(np.sort(x1), np.sort(y1), **kwargs) $\endgroup$ – Igor Rivin Apr 18 '19 at 17:51
  • 1
    $\begingroup$ @IgorRivin excellent improvement! I updated the answer with a modification of your idea and a couple other features. $\endgroup$ – Artem Mavrin Apr 18 '19 at 18:31
  • $\begingroup$ Very nice! Maybe this should be incorporated into the "canon".... $\endgroup$ – Igor Rivin Apr 20 '19 at 18:43
0
$\begingroup$

It might increase your chances of getting a useful answer if you told us your reason for wanting to make such a plot. It will probably be easiest to make a useful plot if your two samples are of the same size.

One possibility is to sort the data and plot one sorted sample against the other. (I'm using R, but I guess you can make sense of the simple code.)

set.seed(1234)
x = rnorm(1000, 100, 10);  y = rexp(1000, 1/100)
plot(sort(x), sort(y), pch=20)

enter image description here

If you have a reasonable guess what the two distributions are (including parameters), you might try transforming each sample by the CDF of its distribution. If your guess is good, you should get (nearly) a straight line.

x1 = pnorm(sort(x), 100, 10)
y1 = pexp(sort(y), 1/100)
plot(x1, y1, pch=20)

enter image description here

If you're trying to see whether the two samples may have come from the same distribution (and in R, the sample sizes are equal and less than 5000), you might use a Kolmogorov-Smirnov Test. (For my x and y, the answer is No.)

ks.test(x, y)

        Two-sample Kolmogorov-Smirnov test

data:  x and y
D = 0.521, p-value < 2.2e-16
alternative hypothesis: two-sided
$\endgroup$
  • $\begingroup$ "Why" is complicated, but in my case I am certain that the distributions are not the same. I would like to know that they are qualitatively of the same "kind". For example, if I look at the eigenvalue spacings of a random hermitian matrix and a random symmetric matrix, the distribution of the spacings will NOT be the same, but is kind of similar, and the q-q plot may capture this. $\endgroup$ – Igor Rivin Apr 17 '19 at 21:12

Not the answer you're looking for? Browse other questions tagged or ask your own question.