Source code for probscale.viz

import copy

import numpy
from matplotlib import pyplot

from .probscale import _minimal_norm
from . import validate
from . import algo

[docs]def probplot(data, ax=None, plottype='prob', dist=None, probax='x', problabel=None, datascale='linear', datalabel=None, bestfit=False, return_best_fit_results=False, estimate_ci=False, ci_kws=None, pp_kws=None, scatter_kws=None, line_kws=None, **fgkwargs): """ Probability, percentile, and quantile plots. Parameters ---------- data : array-like 1-dimensional data to be plotted ax : matplotlib axes, optional The Axes on which to plot. If one is not provided, a new Axes will be created. plottype : string (default = 'prob') Type of plot to be created. Options are: - 'prob': probabilty plot - 'pp': percentile plot - 'qq': quantile plot dist : scipy distribution, optional A distribtion to compute the scale's tick positions. If not specified, a standard normal distribution will be used. probax : string, optional (default = 'x') The axis ('x' or 'y') that will serve as the probability (or quantile) axis. problabel, datalabel : string, optional Axis labels for the probability/quantile and data axes respectively. datascale : string, optional (default = 'log') Scale for the other axis that is not bestfit : bool, optional (default is False) Specifies whether a best-fit line should be added to the plot. return_best_fit_results : bool (default is False) If True a dictionary of results of is returned along with the figure. estimate_ci : bool, optional (False) Estimate and draw a confidence band around the best-fit line using a percentile bootstrap. ci_kws : dict, optional Dictionary of keyword arguments passed directly to ``viz.fit_line`` when computing the best-fit line. pp_kws : dict, optional Dictionary of keyword arguments passed directly to ``viz.plot_pos`` when computing the plotting positions. scatter_kws, line_kws : dict, optional Dictionary of keyword arguments passed directly to ``ax.plot`` when drawing the scatter points and best-fit line, respectively. Other Parameters ---------------- color : string, optional A directly-specified matplotlib color argument for both the data series and the best-fit line if drawn. This argument is made available for compatibility for the seaborn package and is not recommended for general use. Instead colors should be specified within ``scatter_kws`` and ``line_kws``. .. note:: Users should not specify this parameter. It is inteded to only be used by seaborn when operating within a ``FacetGrid``. label : string, optional A directly-specified legend label for the data series. This argument is made available for compatibility for the seaborn package and is not recommended for general use. Instead the data series label should be specified within ``scatter_kws``. .. note:: Users should not specify this parameter. It is inteded to only be used by seaborn when operating within a ``FacetGrid``. Returns ------- fig : matplotlib.Figure The figure on which the plot was drawn. result : dict of linear fit results, optional Keys are: - q : array of quantiles - x, y : arrays of data passed to function - xhat, yhat : arrays of modeled data plotted in best-fit line - res : array of coeffcients of the best-fit line. See also -------- viz.plot_pos viz.fit_line numpy.polyfit scipy.stats.probplot scipy.stats.mstats.plotting_positions Examples -------- Probability plot with the probabilities on the y-axis .. plot:: :context: close-figs >>> import numpy; numpy.random.seed(0) >>> from matplotlib import pyplot >>> from scipy import stats >>> from probscale.viz import probplot >>> data = numpy.random.normal(loc=5, scale=1.25, size=37) >>> fig = probplot(data, plottype='prob', probax='y', ... problabel='Non-exceedance probability', ... datalabel='Observed values', bestfit=True, ... line_kws=dict(linestyle='--', linewidth=2), ... scatter_kws=dict(marker='o', alpha=0.5)) Quantile plot with the quantiles on the x-axis .. plot:: :context: close-figs >>> fig = probplot(data, plottype='qq', probax='x', ... problabel='Theoretical Quantiles', ... datalabel='Observed values', bestfit=True, ... line_kws=dict(linestyle='-', linewidth=2), ... scatter_kws=dict(marker='s', alpha=0.5)) """ if dist is None: dist = _minimal_norm # check input values fig, ax = validate.axes_object(ax) probax = validate.axis_name(probax, 'probability axis') problabel = validate.axis_label(problabel) datalabel = validate.axis_label(datalabel) # default values for symbology options scatter_kws = validate.other_options(scatter_kws) line_kws = validate.other_options(line_kws) pp_kws = validate.other_options(pp_kws) # check plottype plottype = validate.axis_type(plottype) ## !-- kwarg that only seaborn should use --! ## _color = fgkwargs.get('color', None) if _color is not None: scatter_kws['color'] = _color line_kws['color'] = _color ## !-- kwarg that only seaborn should use --! ## _label = fgkwargs.get('label', None) if _label is not None: scatter_kws['label'] = _label # compute the plotting positions and sort the data probs, datavals = plot_pos(data, **pp_kws) qntls = dist.ppf(probs) # determine how the probability values should be expressed if plottype == 'qq': probvals = qntls else: probvals = probs * 100 # set up x, y, Axes for probabilities on the x if probax == 'x': x, y = probvals, datavals ax.set_xlabel(problabel) ax.set_ylabel(datalabel) if plottype == 'prob': ax.set_xscale('prob', dist=dist) fitprobs = 'x' else: fitprobs = None if plottype == 'pp': ax.set_xlim(left=0, right=100) ax.set_yscale(datascale) fitlogs = 'y' if datascale == 'log' else None # setup x, y, Axes for probabilities on the y elif probax == 'y': y, x = probvals, datavals ax.set_xlabel(datalabel) ax.set_ylabel(problabel) if plottype == 'prob': ax.set_yscale('prob', dist=dist) fitprobs = 'y' else: fitprobs = None if plottype == 'pp': ax.set_ylim(bottom=0, top=100) ax.set_xscale(datascale) fitlogs = 'x' if datascale == 'log' else None # finally plot the data linestyle = scatter_kws.pop('linestyle', 'none') marker = scatter_kws.pop('marker', 'o') ax.plot(x, y, linestyle=linestyle, marker=marker, **scatter_kws) # maybe do a best-fit and plot if bestfit: xhat, yhat, modelres = fit_line(x, y, xhat=sorted(x), dist=dist, fitprobs=fitprobs, fitlogs=fitlogs, estimate_ci=estimate_ci) ax.plot(xhat, yhat, **line_kws) if estimate_ci: # for alpha, use half of existing or 0.5 * 0.5 = 0.25 # for zorder, use 1 less than existing or 1 - 1 = 0 opts = { 'facecolor': line_kws.get('color', 'k'), 'edgecolor': line_kws.get('color', 'None'), 'alpha': line_kws.get('alpha', 0.5) * 0.5, 'zorder': line_kws.get('zorder', 1) - 1, 'label': '95% conf. interval' } ax.fill_between(xhat, y1=modelres['yhat_hi'], y2=modelres['yhat_lo'], **opts) else: xhat, yhat, modelres = (None, None, None) # set the probability axes limits if plottype == 'prob': _set_prob_limits(ax, probax, len(probs)) # return the figure and maybe results of the best-fit if return_best_fit_results: results = dict(q=qntls, x=x, y=y, xhat=xhat, yhat=yhat, res=modelres) return fig, results else: return fig
[docs]def plot_pos(data, postype=None, alpha=None, beta=None): """ Compute the plotting positions for a dataset. Heavily borrows from ``scipy.stats.mstats.plotting_positions``. A plottiting position is defined as: ``(i-alpha)/(n+1-alpha-beta)`` where: - ``i`` is the rank order - ``n`` is the size of the dataset - ``alpha`` and ``beta`` are parameters used to adjust the positions. The values of ``alpha`` and ``beta`` can be explicitly set. Typical values can also be access via the ``postype`` parameter. Available ``postype`` values (alpha, beta) are: "type 4" (alpha=0, beta=1) Linear interpolation of the empirical CDF. "type 5" or "hazen" (alpha=0.5, beta=0.5) Piecewise linear interpolation. "type 6" or "weibull" (alpha=0, beta=0) Weibull plotting positions. Unbiased exceedance probability for all distributions. Recommended for hydrologic applications. "type 7" (alpha=1, beta=1) The default values in R. Not recommended with probability scales as the min and max data points get plotting positions of 0 and 1, respectively, and therefore cannot be shown. "type 8" (alpha=1/3, beta=1/3) Approximately median-unbiased. "type 9" or "blom" (alpha=0.375, beta=0.375) Approximately unbiased positions if the data are normally distributed. "median" (alpha=0.3175, beta=0.3175) Median exceedance probabilities for all distributions (used in ``scipy.stats.probplot``). "apl" or "pwm" (alpha=0.35, beta=0.35) Used with probability-weighted moments. "cunnane" (alpha=0.4, beta=0.4) Nearly unbiased quantiles for normally distributed data. This is the default value. "gringorten" (alpha=0.44, beta=0.44) Used for Gumble distributions. Parameters ---------- data : array-like The values whose plotting positions need to be computed. postype : string, optional (default: "cunnane") alpha, beta : float, optional Custom plotting position parameters is the options available through the `postype` parameter are insufficient. Returns ------- plot_pos : numpy.array The computed plotting positions, sorted. data_sorted : numpy.array The original data values, sorted. References ---------- """ pos_params = { 'type 4': (0, 1), 'type 5': (0.5, 0.5), 'type 6': (0, 0), 'type 7': (1, 1), 'type 8': (1/3., 1/3.), 'type 9': (0.375, 0.375), 'weibull': (0, 0), 'median': (0.3175, 0.3175), 'apl': (0.35, 0.35), 'pwm': (0.35, 0.35), 'blom': (0.375, 0.375), 'hazen': (0.5, 0.5), 'cunnane': (0.4, 0.4), 'gringorten': (0.44, 0.44), # Gumble } postype = 'cunnane' if postype is None else postype if alpha is None and beta is None: alpha, beta = pos_params[postype.lower()] data = numpy.asarray(data, dtype=float).flatten() n = data.shape[0] pos = numpy.empty_like(data) pos[n:] = 0 sorted_index = data.argsort() pos[sorted_index[:n]] = (numpy.arange(1, n+1) - alpha) / (n + 1.0 - alpha - beta) return pos[sorted_index], data[sorted_index]
def _set_prob_limits(ax, probax, N): """ Sets the limits of a probabilty axis based the number of point. Parameters ---------- ax : matplotlib Axes The Axes object that will be modified. N : int Maximum number of points for the series plotted on the Axes. which : string The axis whose ticklabels will be rotated. Valid values are 'x', 'y', or 'both'. Returns ------- None """ fig, ax = validate.axes_object(ax) which = validate.axis_name(probax, 'probability axis') if N <= 5: minval = 10 elif N <= 10: minval = 5 else: minval = 10 ** (-1 * numpy.ceil(numpy.log10(N) - 2)) if which in ['x', 'both']: ax.set_xlim(left=minval, right=100-minval) elif which in ['y', 'both']: ax.set_ylim(bottom=minval, top=100-minval)
[docs]def fit_line(x, y, xhat=None, fitprobs=None, fitlogs=None, dist=None, estimate_ci=False, niter=10000, alpha=0.05): """ Fits a line to x-y data in various forms (linear, log, prob scales). Parameters ---------- x, y : array-like Independent and dependent data, respectively. xhat : array-like, optional The values at which ``yhat`` should should be estimated. If not provided, falls back to the sorted values of ``x``. fitprobs, fitlogs : str, optional. Defines how data should be transformed. Valid values are 'x', 'y', or 'both'. If using ``fitprobs``, variables should be expressed as a percentage, i.e., for a probablility transform, data will be transformed with ``lambda x: dist.ppf(x / 100.)``. For a log transform, ``lambda x: numpy.log(x)``. Take care to not pass the same value to both ``fitlogs`` and ``figprobs`` as both transforms will be applied. dist : distribution, optional A fully-spec'd scipy.stats distribution-like object such that ``dist.ppf`` and ``dist.cdf`` can be called. If not provided, defaults to a minimal implementation of ``scipt.stats.norm``. estimate_ci : bool, optional (False) Estimate and draw a confidence band around the best-fit line using a percentile bootstrap. niter : int, optional (default = 10000) Number of bootstrap iterations if ``estimate_ci`` is provided. alpha : float, optional (default = 0.05) The confidence level of the bootstrap estimate. Returns ------- xhat, yhat : numpy arrays Linear model estimates of ``x`` and ``y``. results : dict Dictionary of linear fit results. Keys include: - slope - intersept - yhat_lo (lower confidence interval of the estimated y-vals) - yhat_hi (upper confidence interval of the estimated y-vals) """ fitprobs = validate.fit_argument(fitprobs, "fitprobs") fitlogs = validate.fit_argument(fitlogs, "fitlogs") # maybe set xhat to default values if xhat is None: xhat = copy.copy(x) # maybe set dist to default value if dist is None: dist = _minimal_norm # maybe compute ppf of x if fitprobs in ['x', 'both']: x = dist.ppf(x / 100.) xhat = dist.ppf(numpy.array(xhat)/100.) # maybe compute ppf of y if fitprobs in ['y', 'both']: y = dist.ppf(y / 100.) # maybe compute log of x if fitlogs in ['x', 'both']: x = numpy.log(x) # maybe compute log of y if fitlogs in ['y', 'both']: y = numpy.log(y) yhat, results = algo._fit_simple(x, y, xhat, fitlogs=fitlogs) if estimate_ci: yhat_lo, yhat_hi = algo._bs_fit(x, y, xhat, fitlogs=fitlogs, niter=niter, alpha=alpha) else: yhat_lo, yhat_hi = None, None # maybe undo the ppf transform if fitprobs in ['y', 'both']: yhat = 100. * dist.cdf(yhat) if yhat_lo is not None: yhat_lo = 100. * dist.cdf(yhat_lo) yhat_hi = 100. * dist.cdf(yhat_hi) # maybe undo ppf transform if fitprobs in ['x', 'both']: xhat = 100. * dist.cdf(xhat) results['yhat_lo'] = yhat_lo results['yhat_hi'] = yhat_hi return xhat, yhat, results