knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(pdqr) set.seed(104)
Concept of summary functions is to take one or more pdqr-function(s) and return a summary value (which shouldn't necessarily be a number). Argument method
is used to choose function-specific algorithm of computation.
Note that some summary functions can accumulate pdqr approximation error (like summ_moment()
for example). For better precision increase number intervals for piecewise-linear density using either n
argument for density()
in new_*()
or n_grid
argument in as_*()
.
We will use the following distributions throughout this vignette:
my_beta <- as_d(dbeta, shape1 = 2, shape2 = 5) my_norm <- as_d(dnorm, mean = 0.5) my_beta_mix <- form_mix(list(my_beta, my_beta + 1))
Although they both are continuous, discrete distributions are also fully supported.
# Usage of `summ_center()` summ_center(my_beta, method = "mean") summ_center(my_beta, method = "median") summ_center(my_beta, method = "mode") # Usage of wrappers summ_mean(my_beta) summ_median(my_beta) summ_mode(my_beta) # `summ_mode()` can compute local modes instead of default global summ_mode(my_beta_mix, method = "local")
# Usage of `summ_spread()` summ_spread(my_beta, method = "sd") summ_spread(my_beta, method = "var") summ_spread(my_beta, method = "iqr") summ_spread(my_beta, method = "mad") summ_spread(my_beta, method = "range") # Usage of wrappers summ_sd(my_beta) summ_var(my_beta) summ_iqr(my_beta) summ_mad(my_beta) summ_range(my_beta)
summ_moment()
has extra arguments for controlling the nature of moment (which can be combined):
summ_moment(my_beta, order = 3) summ_moment(my_beta, order = 3, central = TRUE) summ_moment(my_beta, order = 3, standard = TRUE) summ_moment(my_beta, order = 3, absolute = TRUE)
There are wrappers for most common moments: skewness and kurtosis:
summ_skewness(my_beta) # This by default computes excess kurtosis summ_kurtosis(my_beta) # Use `excess = FALSE` to compute non-excess kurtotsis summ_kurtosis(my_beta, excess = FALSE)
summ_quantile(f, probs)
is essentially a more strict version of as_q(f)(probs)
:
summ_quantile(my_beta, probs = seq(0, 1, by = 0.25))
summ_entropy()
computes differential entropy (which can be negative) for "continuous" type pdqr-functions, and information entropy for "discrete":
summ_entropy(my_beta) summ_entropy(new_d(1:10, type = "discrete"))
summ_entropy2()
computes entropy based summary of relation between a pair of distributions. There are two methods: default "relative" (for relative entropy which is Kullback-Leibler divergence) and "cross" (for cross-entropy). It handles different supports by using clip
(default exp(-20)
) value instead of 0 during log()
computation. Order of input does matter: summ_entropy2()
uses support of the first pdqr-function as integration/summation reference.
summ_entropy2(my_beta, my_norm) summ_entropy2(my_norm, my_beta) summ_entropy2(my_norm, my_beta, clip = exp(-10)) summ_entropy2(my_beta, my_norm, method = "cross")
Distributions can be summarized with regions: union of closed intervals. Region is represented as data frame with rows representing intervals and two columns "left" and "right" with left and right interval edges respectively.
summ_interval()
summarizes input pdqr-function with single interval based on the desired coverage level supplied in argument level
. It has three methods:
level
that has minimum width.0.5*(1-level)
and 1 - 0.5*(1-level)
quantiles.level
's critical value (computed from normal distribution). Corresponds to classical confidence interval of sample based on assumption of normality.summ_interval(my_beta, level = 0.9, method = "minwidth") summ_interval(my_beta, level = 0.9, method = "percentile") summ_interval(my_beta, level = 0.9, method = "sigma")
summ_hdr()
computes highest density region (HDR) of a distribution: set of intervals with the lowest total width among all sets with total probability not less than an input level
. With unimodal distribution it is essentially the same as summ_interval()
with "minwidth" method.
# Unimodal distribution summ_hdr(my_beta, level = 0.9) # Multimodal distribution summ_hdr(my_beta_mix, level = 0.9) # Compare this to single interval of minimum width summ_interval(my_beta_mix, level = 0.9, method = "minwidth")
There is a region_*()
family of functions which helps working with them:
beta_mix_hdr <- summ_hdr(my_beta_mix, level = 0.9) beta_mix_interval <- summ_interval(my_beta_mix, level = 0.9) # Test if points are inside region region_is_in(beta_mix_hdr, x = seq(0, 2, by = 0.5)) # Compute total probability of a region region_prob(beta_mix_hdr, f = my_beta_mix) # Pdqr-function doesn't need to be the same as used for computing region region_prob(beta_mix_hdr, f = my_norm) # Compute height of region: minimum value of d-function inside region region_height(beta_mix_hdr, f = my_beta_mix) # Compute width of region: sum of interval widths region_width(beta_mix_hdr) # Compare widths with single interval region_width(beta_mix_interval) # Draw region on existing plot plot(my_beta_mix, main = "90% highest density region") region_draw(beta_mix_hdr)
Function summ_distance()
takes two pdqr-functions and returns a distance between two distributions they represent. Many methods of computation are available. This might be useful for doing comparison statistical inference.
# Kolmogorov-Smirnov distance summ_distance(my_beta, my_norm, method = "KS") # Total variation distance summ_distance(my_beta, my_norm, method = "totvar") # Probability of one distribution being bigger than other, normalized to [0;1] summ_distance(my_beta, my_norm, method = "compare") # Wassertein distance: "average path density point should travel while # transforming from one into another" summ_distance(my_beta, my_norm, method = "wass") # Cramer distance: integral of squared difference of p-functions summ_distance(my_beta, my_norm, method = "cramer") # "Align" distance: path length for which one of distribution should be "moved" # towards the other so that they become "aligned" (probability of one being # greater than the other is 0.5) summ_distance(my_beta, my_norm, method = "align") # "Entropy" distance: `KL(f, g) + KL(g, f)`, where `KL()` is Kullback-Leibler # divergence. Usually should be used for distributions with same support, but # works even if they are different (with big numerical penalty). summ_distance(my_beta, my_norm, method = "entropy")
Function summ_separation()
computes a threshold that optimally separates distributions represented by pair of input pdqr-functions. In other words, summ_separation()
solves a binary classification problem with one-dimensional linear classifier: values not more than some threshold are classified as one class, and more than threshold - as another. Order of input functions doesn't matter.
summ_separation(my_beta, my_norm, method = "KS") summ_separation(my_beta, my_norm, method = "F1")
Functions summ_classmetric()
and summ_classmetric_df()
compute metric(s) of classification setup, similar to one used in summ_separation()
. Here classifier threshold should be supplied and order of input matters. Classification is assumed to be done as follows: any x value not more than threshold value is classified as "negative"; if more - "positive". Classification metrics are computed based on two pdqr-functions: f
, which represents the distribution of values which should be classified as "negative" ("true negative"), and g
- the same for "positive" ("true positive").
# Many threshold values can be supplied thres_vec <- seq(0, 1, by = 0.2) summ_classmetric(f = my_beta, g = my_norm, threshold = thres_vec, method = "F1") # In `summ_classmetric_df()` many methods can be supplied summ_classmetric_df( f = my_beta, g = my_norm, threshold = thres_vec, method = c("GM", "F1", "MCC") )
With summ_roc()
and summ_rocauc()
one can compute data frame of ROC curve points and ROC AUC value respectively. There is also a roc_plot()
function for predefined plotting of ROC curve.
my_roc <- summ_roc(my_beta, my_norm) head(my_roc) summ_rocauc(my_beta, my_norm) roc_plot(my_roc)
'pdqr' has functions that can order set of distributions. They are summ_order()
, summ_sort()
, and summ_rank()
, which are analogues of order()
, sort()
, and rank()
respectively. They take a list of pdqr-functions as input, establish their ordering based on specified method, and return the desired output.
There are two sets of methods:
f
is greater than g
if and only if P(f >= g) > 0.5
, or in 'pdqr' code summ_prob_true(f >= g) > 0.5
. This method orders input based on this relation and order()
function. Notes:order()
.summ_center()
: ordering of distributions is defined as ordering of corresponding measures of distribution's center.# Here the only clear "correct" ordering is that `a <= b`. f_list <- list(a = my_beta, b = my_beta + 1, c = my_norm) # Returns an integer vector representing a permutation which rearranges f_list # in desired order summ_order(f_list, method = "compare") # In this particular case of `f_list` all orderings agree with each other, but # generally this is not the case: for any pair of methods there is a case # when they disagree with each other summ_order(f_list, method = "mean") summ_order(f_list, method = "median") summ_order(f_list, method = "mode") # Use `decreasing = TRUE` to sort decreasingly summ_order(f_list, method = "compare", decreasing = TRUE) # Sort list summ_sort(f_list) summ_sort(f_list, decreasing = TRUE) # Rank elements: 1 indicates "the smallest", `length(f_list)` - "the biggest" summ_rank(f_list)
Functions summ_prob_true()
and summ_prob_false()
should be used to extract probabilities from boolean pdqr-functions: outputs of comparing basic operators (like >=
, ==
, etc.):
summ_prob_true(my_beta >= my_norm) summ_prob_false(my_beta >= 2*my_norm)
summ_pval()
computes p-value(s) of observed statistic(s) based on the distribution. You can compute left, right, or two-sided p-values with methods "left", "right", and "both" respectively. By default multiple input values are adjusted for multiple comparisons (using stats::p.adjust()):
# By default two-sided p-value is computed summ_pval(my_beta, obs = 0.7) summ_pval(my_beta, obs = 0.7, method = "left") summ_pval(my_beta, obs = 0.7, method = "right") # Multiple values are adjusted with `p.adjust()` with "holm" method by default obs_vec <- seq(0, 1, by = 0.1) summ_pval(my_beta, obs = obs_vec) # Use `adjust = "none"` to not adjust summ_pval(my_beta, obs = obs_vec, adjust = "none")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.