Source code for dendropy.calculate.statistics

#! /usr/bin/env python

##############################################################################
##  DendroPy Phylogenetic Computing Library.
##
##  Copyright 2010-2015 Jeet Sukumaran and Mark T. Holder.
##  All rights reserved.
##
##  See "LICENSE.rst" for terms and conditions of usage.
##
##  If you use this work or any portion thereof in published work,
##  please cite it as:
##
##     Sukumaran, J. and M. T. Holder. 2010. DendroPy: a Python library
##     for phylogenetic computing. Bioinformatics 26: 1569-1571.
##
##############################################################################

"""
Functions to calculate some general statistics.
"""

import math
from dendropy.calculate import probability
from operator import itemgetter

def _mean_and_variance_pop_n(values):
    n = 0
    s = 0.0
    ss = 0.0
    for v in values:
        n += 1
        s += v
        ss += v*v
    if n == 0:
        raise IndexError("values in mean_and_variance cannot be empty")
    mean = float(s)/n
    var = (ss - mean*s)/n
    return mean, var, n

[docs]def mean_and_population_variance(values): """Returns the mean and population variance while only passing over the elements in values once.""" return _mean_and_variance_pop_n(values)[:2]
[docs]def mean_and_sample_variance(values): """Returns the mean and sample variance while only passing over the elements in values once.""" mean, pop_var, n = _mean_and_variance_pop_n(values) if n == 1: return mean, float('inf') samp_var = n*pop_var/(n-1) return mean, samp_var
[docs]def mode(values, bin_size=0.1): """ Returns the mode of a set of values. """ bins = {} for v in values: if bin_size is not None: idx = int(round(float(v)/bin_size)) else: idx = v if idx in bins: bins[idx] += 1 else: bins[idx] = 1 sorted_bins = sorted(bins.items(), key=itemgetter(1), reverse=True) max_count = sorted_bins[0][1] results = [(sorted_bins[i][0] * bin_size) for i in xrange(len(sorted_bins)) if sorted_bins[i][1] >= max_count] return results
[docs]def median(pool): """ Returns median of sample. From: http://wiki.python.org/moin/SimplePrograms """ copy = sorted(pool) size = len(copy) if size % 2 == 1: idx = int((size - 1) / 2) return copy[idx] else: idx1 = int(size/2) - 1 idx2 = int(size/2) return (copy[idx1] + copy[idx2]) / 2
[docs]def empirical_hpd(values, conf=0.05): """ Assuming a **unimodal** distribution, returns the 0.95 highest posterior density (HPD) interval for a set of samples from a posterior distribution. Adapted from ``emp.hpd`` in the "TeachingDemos" R package (Copyright Greg Snow; licensed under the Artistic License). """ conf = min([conf, 1.0 - conf]) n = len(values) nn = int(round(n * conf)) x = sorted(values) xx = [] if nn == 0: raise ValueError("Sample size too small: %s" % len(values)) for i in range(nn): Z1 = x[n-(nn-i)] Z2 = x[i] #print "==>", Z1, Z2 xx.append(Z1 - Z2) m = min(xx) nnn = xx.index(m) #print "n=", n #print "nn=", nn #print "xx=", xx #print "m=", m #print "nnn=", n return (x[nnn], x[n-nn+nnn])
[docs]def empirical_cdf(values, v): """ Returns the proportion of values in ``values`` <= ``v``. """ count = 0.0 for idx, v0 in enumerate(values): if v0 < v: count += 1 return count / len(values)
def quantile(values, q): """ Returns q-th quantile. """ values = sorted(values) size = len(values) idx = int(round(size * q)) - 1 if idx == 0: raise ValueError("Sample size too small: %s" % len(values)) return values[idx] # http://adorio-research.org/wordpress/?p=125 # File quantile.py # Desc computes sample quantiles # Author Ernesto P. Adorio, PhD. # UPDEPP (U.P. at Clarkfield) # Version 0.0.1 August 7. 2009 def quantile(x, q, qtype = 7, issorted = False): from math import modf, floor """ Args: x - input data q - quantile qtype - algorithm issorted- True if x already sorted. Compute quantiles from input array x given q.For median, specify q=0.5. References: http://reference.wolfram.com/mathematica/ref/Quantile.html http://wiki.r-project.org/rwiki/doku.php?id=rdoc:stats:quantile Author: Ernesto P.Adorio Ph.D. UP Extension Program in Pampanga, Clark Field. """ if not issorted: y = sorted(x) else: y = x if not (1 <= qtype <= 9): return None # error! # Parameters for the Hyndman and Fan algorithm abcd = [(0, 0, 1, 0), # inverse empirical distrib.function., R type 1 (0.5, 0, 1, 0), # similar to type 1, averaged, R type 2 (0.5, 0, 0, 0), # nearest order statistic,(SAS) R type 3 (0, 0, 0, 1), # California linear interpolation, R type 4 (0.5, 0, 0, 1), # hydrologists method, R type 5 (0, 1, 0, 1), # mean-based estimate(Weibull method), (SPSS,Minitab), type 6 (1, -1, 0, 1), # mode-based method,(S, S-Plus), R type 7 (1.0/3, 1.0/3, 0, 1), # median-unbiased , R type 8 (3/8.0, 0.25, 0, 1) # normal-unbiased, R type 9. ] a, b, c, d = abcd[qtype-1] n = len(x) g, j = modf( a + (n+b) * q -1) if j < 0: return y[0] elif j >= n: return y[n-1] # oct. 8, 2010 y[n]???!! uncaught off by 1 error!!! j = int(floor(j)) if g == 0: return y[j] else: return y[j] + (y[j+1]- y[j])* (c + d * g)
[docs]def quantile_5_95(values): """ Returns 5% and 95% quantiles. """ values = sorted(values) size = len(values) idx5 = int(round(size * 0.05)) - 1 idx95 = int(round(size * 0.95)) - 1 if idx5 == 0: raise ValueError("Sample size too small: %s" % len(values)) return values[idx5], values[idx95]
[docs]def variance_covariance(data, population_variance=False): """ Returns the Variance-Covariance matrix for ``data``. From: http://www.python-forum.org/pythonforum/viewtopic.php?f=3&t=17441 """ N = len(data) # number of vectors D = len(data[0]) # dimensions per vector if population_variance: denom = N else: denom = N-1.0 means = [0.0 for i in range(D)] # intialize 1xD mean vector for i in range(N): for j in range(D): means[j] += data[i][j] means = [i/N for i in means] # print "Means:"," ".join(map(str,means)),"\n" covar = [[0.0 for i in range(D)] for j in range(D)] # initialize DxD covariance matrix for i in range(D): for j in range(i+1): # covariance symmetric, only do lower triangle of matrix sum = 0.0 for k in range(N): sum += data[k][i]*data[k][j] covar[i][j] = sum/denom - means[i]*means[j]*N/denom for j in range(D): for k in range(j+1): covar[k][j] = covar[j][k] # print "covariance:" # for i in range(D): # print " ".join(map(str,covar[i])) # print "" return covar
[docs]def rank(value_to_be_ranked, value_providing_rank): """ Returns the rank of ``value_to_be_ranked`` in set of values, ``values``. Works even if ``values`` is a non-orderable collection (e.g., a set). A binary search would be an optimized way of doing this if we can constrain ``values`` to be an ordered collection. """ num_lesser = [v for v in value_providing_rank if v < value_to_be_ranked] return len(num_lesser)
[docs]class FishersExactTest(object): """ Given a 2x2 table: +---+---+ | a | b | +---+---+ | c | d | +---+---+ represented by a list of lists:: [[a,b],[c,d]] this calculates the sum of the probability of this table and all others more extreme under the null hypothesis that there is no association between the categories represented by the vertical and horizontal axes. """
[docs] def probability_of_table(table): """ Given a 2x2 table: +---+---+ | a | b | +---+---+ | c | d | +---+---+ represented by a list of lists:: [[a,b],[c,d]] this returns the probability of this table under the null hypothesis of no association between rows and columns, which was shown by Fisher to be a hypergeometric distribution: p = ( choose(a+b, a) * choose(c+d, c) ) / choose(a+b+c+d, a+c) """ a = table[0][0] b = table[0][1] c = table[1][0] d = table[1][1] return probability.hypergeometric_pmf(a, a+b, c+d, a+c)
probability_of_table = staticmethod(probability_of_table) def __init__(self, table): self.table = table self.flat_table = [table[0][0], table[0][1], table[1][0], table[1][1]] self.min_value = min(self.flat_table) self.max_value = max(self.flat_table) def _rotate_cw(self, table): """ Returns a copy of table such that all the values are rotated clockwise once. """ return [ [ table[1][0], table[0][0] ], [table[1][1], table[0][1] ] ] def _min_rotation(self): """ Returns copy of self.table such that the smallest value is in the first (upper left) cell. """ table = [list(self.table[0]), list(self.table[1])] while table[0][0] != self.min_value: table = self._rotate_cw(table) return table def _max_rotation(self): """ Returns copy of self.table such that the largest value is in the first (upper left) cell. """ table = [list(self.table[0]), list(self.table[1])] while table[0][0] != self.max_value: table = self._rotate_cw(table) return table def _sum_left_tail(self): """ Returns the sum of probabilities of tables that are *more* extreme than the current table. """ # left_tail_tables = self._get_left_tail_tables() # p_vals = [ self.probability_of_table(t) for t in left_tail_tables ] p_vals = self._get_left_tail_probs() return sum(p_vals) def _sum_right_tail(self): """ Returns the sum of probabilities of tables that are *less* extreme than the current table. """ # right_tail_tables = self._get_right_tail_tables() # p_vals = [ self.probability_of_table(t) for t in right_tail_tables ] p_vals = self._get_right_tail_probs() return sum(p_vals) def _get_left_tail_probs(self): """ Returns list of probabilities of all tables *more* extreme than the current table. """ table = self._min_rotation() row_totals = [sum(table[0]), sum(table[1])] col_totals = [table[0][0] + table[1][0], table[0][1] + table[1][1]] p_vals = [] while True: table[0][0] -= 1 if table[0][0] < 0: break table[0][1] = row_totals[0] - table[0][0] table[1][0] = col_totals[0] - table[0][0] table[1][1] = row_totals[1] - table[1][0] p_vals.append(self.probability_of_table(table)) return p_vals def _get_right_tail_probs(self): """ Returns list of probabilities of all tables *less* extreme than the current table. """ table = self._min_rotation() row_totals = [sum(table[0]), sum(table[1])] col_totals = [table[0][0] + table[1][0], table[0][1] + table[1][1]] p_vals = [] while True: table[0][0] += 1 table[0][1] = row_totals[0] - table[0][0] if table[0][1] < 0: break table[1][0] = col_totals[0] - table[0][0] if table[1][0] < 0: break table[1][1] = row_totals[1] - table[1][0] if table[1][1] < 0: break p_vals.append(self.probability_of_table(table)) return p_vals def _get_left_tail_tables(self): """ Returns all tables that are *more* extreme than the current table. """ table = self._min_rotation() row_totals = [sum(table[0]), sum(table[1])] col_totals = [table[0][0] + table[1][0], table[0][1] + table[1][1]] left_tail_tables = [] while True: table[0][0] -= 1 if table[0][0] < 0: break table[0][1] = row_totals[0] - table[0][0] table[1][0] = col_totals[0] - table[0][0] table[1][1] = row_totals[1] - table[1][0] left_tail_tables.append([list(table[0]), list(table[1])]) return left_tail_tables def _get_right_tail_tables(self): """ Returns all tables that are *less* extreme than the current table. """ table = self._min_rotation() row_totals = [sum(table[0]), sum(table[1])] col_totals = [table[0][0] + table[1][0], table[0][1] + table[1][1]] right_tail_tables = [] while True: table[0][0] += 1 table[0][1] = row_totals[0] - table[0][0] if table[0][1] < 0: break table[1][0] = col_totals[0] - table[0][0] if table[1][0] < 0: break table[1][1] = row_totals[1] - table[1][0] if table[1][1] < 0: break right_tail_tables.append([list(table[0]), list(table[1])]) return right_tail_tables
[docs] def left_tail_p(self): """ Returns the sum of probabilities of this table and all others more extreme. """ return self.probability_of_table(self.table) + self._sum_left_tail()
[docs] def right_tail_p(self): """ Returns the sum of probabilities of this table and all others more extreme. """ return self.probability_of_table(self.table) + self._sum_right_tail()
[docs] def two_tail_p(self): """ Returns the sum of probabilities of this table and all others more extreme. """ p0 = self.probability_of_table(self.table) all_p_vals = self._get_left_tail_probs() + self._get_right_tail_probs() p_vals = [] for p in all_p_vals: if p <= p0: p_vals.append(p) return sum(p_vals) + p0
[docs]def summarize(values): """ Summarizes a sample of values: - ``range`` : tuple pair representing minimum and maximum values - ``mean`` : mean of sample - ``median`` : median of sample - ``var`` : (sample) variance - ``sd`` : (sample) standard deviation - ``hpd95`` : tuple pair representing 5% and 95% HPD - ``quant_5_95`` : tuple pair representing 5% and 95% quantile """ summary = {} if len(values) == 0: raise ValueError("No values in data") try: summary['range'] = (min(values), max(values)) except (ValueError, OverflowError): summary['range'] = None try: summary['mean'], summary['var'] = mean_and_sample_variance(values) try: #summary['sd'] = math.sqrt(summary['var']) summary['sd'] = summary['var'] ** 0.5 except ValueError: summary['sd'] = 0.0 except OverflowError: summary['sd'] = None except (ValueError, OverflowError, IndexError): summary['mean'], summary['var'], summary['sd'] = None, None, None try: summary['median'] = median(values) except (ValueError, OverflowError): summary['median'] = None try: summary['hpd95'] = empirical_hpd(values, conf=0.95) except (ValueError, OverflowError): summary['hpd95'] = None try: summary['quant_5_95'] = quantile_5_95(values) except (ValueError, OverflowError): summary['quant_5_95'] = None return summary