basic statistics module this module provides functions for calculating statistics of data including averages variance and standard deviation calculating averages function description mean arithmetic mean average of data fmean fast floating point arithmetic mean geometricmean geometric mean of data harmonicmean harmonic mean of data median median middle value of data medianlow low median of data medianhigh high median of data mediangrouped median or 50th percentile of grouped data mode mode most common value of data multimode list of modes most common values of data quantiles divide data into intervals with equal probability calculate the arithmetic mean the average of data mean1 0 2 5 3 25 5 75 2 625 calculate the standard median of discrete data median2 3 4 5 3 5 calculate the median or 50th percentile of data grouped into class intervals centred on the data values provided e g if your data points are rounded to the nearest whole number mediangrouped2 2 3 3 3 4 doctest ellipsis 2 8333333333 this should be interpreted in this way you have two data points in the class interval 1 52 5 three data points in the class interval 2 53 5 and one in the class interval 3 54 5 the median of these data points is 2 8333 calculating variability or spread function description pvariance population variance of data variance sample variance of data pstdev population standard deviation of data stdev sample standard deviation of data calculate the standard deviation of sample data stdev2 5 3 25 5 5 11 25 11 75 doctest ellipsis 4 38961843444 if you have previously calculated the mean you can pass it as the optional second argument to the four spread functions to avoid recalculating it data 1 2 2 4 4 4 5 6 mu meandata pvariancedata mu 2 5 statistics for relations between two inputs function description covariance sample covariance for two variables correlation pearson s correlation coefficient for two variables linearregression intercept and slope for simple linear regression calculate covariance pearson s correlation and simple linear regression for two inputs x 1 2 3 4 5 6 7 8 9 y 1 2 3 1 2 3 1 2 3 covariancex y 0 75 correlationx y doctest ellipsis 0 31622776601 linearregressionx y doctest linearregressionslope0 1 intercept1 5 exceptions a single exception is defined statisticserror is a subclass of valueerror exceptions private utilities sumdata type sum count return a highprecision sum of the given numeric data as a fraction together with the type to be converted to and the count of items examples sum3 2 25 4 5 0 5 0 25 class float fraction19 2 5 some sources of roundoff error will be avoided builtin sum returns zero sum1e50 1 1e50 1000 class float fraction1000 1 3000 fractions and decimals are also supported from fractions import fraction as f sumf2 3 f7 5 f1 4 f5 6 class fractions fraction fraction63 20 4 from decimal import decimal as d data d0 1375 d0 2108 d0 3061 d0 0419 sumdata class decimal decimal fraction6963 10000 4 mixed types are currently treated as an error except that int is allowed the sum will be a nan or inf we can ignore all the finite partials and just look at this special one sum all the partial sums using builtin sum return the exact mean and sum of square deviations of sequence data calculations are done in a single pass allowing the input to be an iterator if given c is used the mean otherwise it is calculated from the data use the c argument with care as it can lead to garbage results the sum will be a nan or inf we can ignore all the finite partials and just look at this special one this formula has poor numeric properties for floats but with fractions it is exact coerce types t and s to a common type or raise typeerror coercion rules are currently an implementation detail see the coercetest test class in teststatistics for details see http bugs python orgissue24068 if the types are the same no need to coerce anything put this first so that the usual case no coercion needed happens as soon as possible mixed int other coerce to the other type if one is a strict subclass of the other coerce to the subclass ints coerce to the other type mixed fraction float coerces to float or float subclass any other combination is disallowed return real number x to exact numerator denominator pair exactratio0 25 1 4 x is expected to be an int fraction decimal or float xxx we should revisit whether using fractions to accumulate exact ratios is the right way to go the integer ratios for binary floats can have numerators or denominators with over 300 decimal digits the problem is more acute with decimal floats where the default decimal context supports a huge range of exponents from emin999999 to emax999999 when expanded with asintegerratio numbers like decimal 3 14e5000 and decimal 3 14e5000 have large numerators or denominators that will slow computation when the integer ratios are accumulated as fractions the size grows to cover the full range from the smallest magnitude to the largest for example fraction3 14e300 fraction3 14e300 has a 616 digit numerator likewise fractiondecimal 3 14e5000 fractiondecimal 3 14e5000 has 10 003 digit numerator this doesn t seem to have been problem in practice but it is a potential pitfall float nan or inf x may be an integral abc convert value to given numeric type t if typevalue is t this covers the cases where t is fraction or where value is a nan or inf decimal or float return value if issubclasst int and value denominator 1 t float try fixme what do we do if this overflows return tvalue except typeerror if issubclasst decimal return tvalue numerator tvalue denominator else raise def failnegvalues errmsg negative value rank order a dataset the lowest value has rank 1 ties are averaged so that equal values receive the same rank data 31 56 31 25 75 18 rankdata 3 5 5 0 3 5 2 0 6 0 1 0 the operation is idempotent rank3 5 5 0 3 5 2 0 6 0 1 0 3 5 5 0 3 5 2 0 6 0 1 0 it is possible to rank the data in reverse order so that the highest value has rank 1 also a keyfunction can extract the field to be ranked goals eagles 45 bears 48 lions 44 rankgoals keyitemgetter1 reversetrue 2 0 1 0 3 0 ranks are conventionally numbered starting from one however setting start to zero allows the ranks to be used as array indices prize gold silver bronze certificate scores 8 1 7 3 9 4 8 3 prizeinti for i in rankscores start0 reversetrue bronze certificate gold silver if this function becomes public at some point more thought needs to be given to the signature a list of ints is plausible when ties is min or max when ties is average either listfloat or listfraction is plausible default handling of ties matches scipy stats mstats spearmanr square root of nm rounded to the nearest integer using roundtoodd reference https www lri frmelquiondoc05imacs171expose pdf a math isqrtn m return a aam n for 53 bit precision floats the bit width used in floatsqrtoffrac is 109 sqrtbitwidth int 2 sys floatinfo mantdig 3 def floatsqrtoffracn int m int float see principle and proof sketch at https bugs python orgmsg407078 square root of nm as a decimal correctly rounded premise for decimal computing nm sqrt can be off by 1 ulp from the correctly rounded result method check the result moving up or down a step if needed if n 0 if not n return decimal 0 0 n m n m root decimaln decimalm sqrt nr dr root asintegerratio plus root nextplus np dp plus asintegerratio test n m root plus 2 2 if 4 n drdp2 m drnp dpnr2 return plus minus root nextminus nm dm minus asintegerratio test n m root minus 2 2 if 4 n drdm2 m drnm dmnr2 return minus return root measures of central tendency averages def meandata t total n sumdata if n 1 raise statisticserror mean requires at least one data point return converttotal n t def fmeandata weightsnone if weights is none try n lendata except typeerror handle iterators that do not define len n 0 def countiterable nonlocal n for n x in enumerateiterable start1 yield x data countdata total fsumdata if not n raise statisticserror fmean requires at least one data point return total n if not isinstanceweights list tuple weights listweights try num sumproddata weights except valueerror raise statisticserror data and weights must be the same length den fsumweights if not den raise statisticserror sum of weights must be nonzero return num den def geometricmeandata try return expfmeanmaplog data except valueerror raise statisticserror geometric mean requires a nonempty dataset containing positive numbers from none def harmonicmeandata weightsnone if iterdata is data data listdata errmsg harmonic mean does not support negative values n lendata if n 1 raise statisticserror harmonicmean requires at least one data point elif n 1 and weights is none x data0 if isinstancex numbers real decimal if x 0 raise statisticserrorerrmsg return x else raise typeerror unsupported type if weights is none weights repeat1 n sumweights n else if iterweights is weights weights listweights if lenweights n raise statisticserror number of weights does not match data size sumweights sumw for w in failnegweights errmsg try data failnegdata errmsg t total count sumw x if w else 0 for w x in zipweights data except zerodivisionerror return 0 if total 0 raise statisticserror weighted sum must be positive return convertsumweights total t fixme investigate ways to calculate medians without sorting quickselect def mediandata data sorteddata n lendata if n 0 raise statisticserrorno median for empty data if n 2 1 return datan 2 else i n 2 return datai 1 datai 2 def medianlowdata data sorteddata n lendata if n 0 raise statisticserrorno median for empty data if n 2 1 return datan 2 else return datan 2 1 def medianhighdata data sorteddata n lendata if n 0 raise statisticserrorno median for empty data return datan 2 def mediangroupeddata interval1 0 data sorteddata n lendata if not n raise statisticserrorno median for empty data find the value at the midpoint remember this corresponds to the midpoint of the class interval x datan 2 using olog n bisection find where all the x values occur in the data all x will lie within datai j i bisectleftdata x j bisectrightdata x loi coerce to floats raising a typeerror if not possible try interval floatinterval x floatx except valueerror raise typeerrorf value cannot be converted to a float interpolate the median using the formula found at https www cuemath comdatamedianofgroupeddata l x interval 2 0 lower limit of the median interval cf i cumulative frequency of the preceding interval f j i number of elements in the median internal return l interval n 2 cf f def modedata pairs counteriterdata mostcommon1 try return pairs00 except indexerror raise statisticserror no mode for empty data from none def multimodedata counts counteriterdata if not counts return maxcount maxcounts values return value for value count in counts items if count maxcount notes on methods for computing quantiles there is no one perfect way to compute quantiles here we offer two methods that serve common needs most other packages surveyed offered at least one or both of these two making them standard in the sense of widelyadopted and reproducible they are also easy to explain easy to compute manually and have straightforward interpretations that aren t surprising the default method is known as r6 percentile exc or expected value of rank order statistics the alternative method is known as r7 percentile inc or mode of rank order statistics for sample data where there is a positive probability for values beyond the range of the data the r6 exclusive method is a reasonable choice consider a random sample of nine values from a population with a uniform distribution from 0 0 to 1 0 the distribution of the third ranked sample point is described by betavariatealpha3 beta7 which has mode0 250 median0 286 and mean0 300 only the latter which corresponds with r6 gives the desired cut point with 30 of the population falling below that value making it comparable to a result from an invcdf function the r6 exclusive method is also idempotent for describing population data where the end points are known to be included in the data the r7 inclusive method is a reasonable choice instead of the mean it uses the mode of the beta distribution for the interior points per hyndman fan one nice property is that the vertices of q7p divide the range into n 1 intervals and exactly 100p of the intervals lie to the left of q7p and 1001 p of the intervals lie to the right of q7p if needed other methods could be added however for now the position is that fewer options make for easier choices and that external packages can be used for anything more advanced def quantilesdata n4 method exclusive if n 1 raise statisticserror n must be at least 1 data sorteddata ld lendata if ld 2 if ld 1 return data n 1 raise statisticserror must have at least one data point if method inclusive m ld 1 result for i in range1 n j delta divmodi m n interpolated dataj n delta dataj 1 delta n result appendinterpolated return result if method exclusive m ld 1 result for i in range1 n j i m n rescale i to mn j 1 if j 1 else ld1 if j ld1 else j clamp to 1 ld1 delta im jn exact integer math interpolated dataj 1 n delta dataj delta n result appendinterpolated return result raise valueerrorf unknown method method r measures of spread see http mathworld wolfram comvariance html http mathworld wolfram comsamplevariance html def variancedata xbarnone t ss c n ssdata xbar if n 2 raise statisticserror variance requires at least two data points return convertss n 1 t def pvariancedata munone t ss c n ssdata mu if n 1 raise statisticserror pvariance requires at least one data point return convertss n t def stdevdata xbarnone t ss c n ssdata xbar if n 2 raise statisticserror stdev requires at least two data points mss ss n 1 if issubclasst decimal return decimalsqrtoffracmss numerator mss denominator return floatsqrtoffracmss numerator mss denominator def pstdevdata munone t ss c n ssdata mu if n 1 raise statisticserror pstdev requires at least one data point mss ss n if issubclasst decimal return decimalsqrtoffracmss numerator mss denominator return floatsqrtoffracmss numerator mss denominator def meanstdevdata handle nans and infs gracefully finite inputs overflowed so scale down and recompute nonzero inputs underflowed so scale up and recompute scale 1 sqrtsys floatinfo min sys floatinfo epsilon improve accuracy with a differential correction https www wolframalpha cominput imaclaurinseriessqrt28h22bx29atx3d0 statistics for relations between two inputs see https en wikipedia orgwikicovariance https en wikipedia orgwikipearsoncorrelationcoefficient https en wikipedia orgwikisimplelinearregression covariance return the sample covariance of two inputs x and y covariance is a measure of the joint variability of two inputs x 1 2 3 4 5 6 7 8 9 y 1 2 3 1 2 3 1 2 3 covariancex y 0 75 z 9 8 7 6 5 4 3 2 1 covariancex z 7 5 covariancez x 7 5 pearson s correlation coefficient return the pearson s correlation coefficient for two inputs pearson s correlation coefficient r takes values between 1 and 1 it measures the strength and direction of a linear relationship x 1 2 3 4 5 6 7 8 9 y 9 8 7 6 5 4 3 2 1 correlationx x 1 0 correlationx y 1 0 if method is ranked computes spearman s rank correlation coefficient for two inputs the data is replaced by ranks ties are averaged so that equal values receive the same rank the resulting coefficient measures the strength of a monotonic relationship spearman s rank correlation coefficient is appropriate for ordinal data or for continuous data that doesn t meet the linear proportion requirement for pearson s correlation coefficient slope and intercept for simple linear regression return the slope and intercept of simple linear regression parameters estimated using ordinary least squares simple linear regression describes relationship between an independent variable x and a dependent variable y in terms of a linear function y slope x intercept noise where slope and intercept are the regression parameters that are estimated and noise represents the variability of the data that was not explained by the linear regression it is equal to the difference between predicted and actual values of the dependent variable the parameters are returned as a named tuple x 1 2 3 4 5 noise normaldist samples5 seed42 y 3 xi 2 noisei for i in range5 linearregressionx y doctest ellipsis linearregressionslope3 17495 intercept1 00925 if proportional is true the independent variable x and the dependent variable y are assumed to be directly proportional the data is fit to a line passing through the origin since the intercept will always be 0 0 the underlying linear function simplifies to y slope x noise y 3 xi noisei for i in range5 linearregressionx y proportionaltrue doctest ellipsis linearregressionslope2 90475 intercept0 0 normal distribution there is no closedform solution to the inverse cdf for the normal distribution so we use a rational approximation instead wichura m j 1988 algorithm as241 the percentage points of the normal distribution applied statistics blackwell publishing 37 3 477484 doi 10 23072347330 jstor 2347330 hash sum 55 8831928806149014439 hash sum 49 33206503301610289036 hash sum 47 52583317549289671629 if available use c implementation https en wikipedia orgwikinormaldistribution https en wikipedia orgwikivarianceproperties inverse cumulative distribution function x px x p finds the value of the random variable such that the probability of the variable being less than or equal to that value equals the given probability this function is also called the percent point function or quantile function divide into n continuous intervals with equal probability returns a list of n 1 cut points separating the intervals set n to 4 for quartiles the default set n to 10 for deciles set n to 100 for percentiles which gives the 99 cuts points that separate the normal distribution in to 100 equal sized groups compute the overlapping coefficient ovl between two normal distributions measures the agreement between two normal probability distributions returns a value between 0 0 and 1 0 giving the overlapping area in the two underlying probability density functions n1 normaldist2 4 1 6 n2 normaldist3 2 2 0 n1 overlapn2 0 8035050657330205 see the overlapping coefficient as a measure of agreement between probability distributions and point estimation of the overlap of two normal densities henry f inman and edwin l bradley jr http dx doi org10 108003610928908830127 compute the standard score x mean stdev describes x in terms of the number of standard deviations above or below the mean of the normal distribution https www statisticshowto comprobabilityandstatisticszscore return the mode of the normal distribution the mode is the value x where which the probability density function pdf takes its maximum value add a constant or another normaldist instance if other is a constant translate mu by the constant leaving sigma unchanged if other is a normaldist add both the means and the variances mathematically this works only if the two distributions are independent or if they are jointly normally distributed subtract a constant or another normaldist instance if other is a constant translate by the constant mu leaving sigma unchanged if other is a normaldist subtract the means and add the variances mathematically this works only if the two distributions are independent or if they are jointly normally distributed multiply both mu and sigma by a constant used for rescaling perhaps to change measurement units sigma is scaled with the absolute value of the constant divide both mu and sigma by a constant used for rescaling perhaps to change measurement units sigma is scaled with the absolute value of the constant exceptions private utilities _sum data type sum count return a high precision sum of the given numeric data as a fraction together with the type to be converted to and the count of items examples _sum 3 2 25 4 5 0 5 0 25 class float fraction 19 2 5 some sources of round off error will be avoided built in sum returns zero _sum 1e50 1 1e50 1000 class float fraction 1000 1 3000 fractions and decimals are also supported from fractions import fraction as f _sum f 2 3 f 7 5 f 1 4 f 5 6 class fractions fraction fraction 63 20 4 from decimal import decimal as d data d 0 1375 d 0 2108 d 0 3061 d 0 0419 _sum data class decimal decimal fraction 6963 10000 4 mixed types are currently treated as an error except that int is allowed the sum will be a nan or inf we can ignore all the finite partials and just look at this special one sum all the partial sums using builtin sum or raise typeerror return the exact mean and sum of square deviations of sequence data calculations are done in a single pass allowing the input to be an iterator if given c is used the mean otherwise it is calculated from the data use the c argument with care as it can lead to garbage results the sum will be a nan or inf we can ignore all the finite partials and just look at this special one this formula has poor numeric properties for floats but with fractions it is exact or raise typeerror likely a decimal coerces to float first coerce types t and s to a common type or raise typeerror coercion rules are currently an implementation detail see the coercetest test class in test_statistics for details see http bugs python org issue24068 if the types are the same no need to coerce anything put this first so that the usual case no coercion needed happens as soon as possible mixed int other coerce to the other type if one is a strict subclass of the other coerce to the subclass ints coerce to the other type mixed fraction float coerces to float or float subclass any other combination is disallowed return real number x to exact numerator denominator pair _exact_ratio 0 25 1 4 x is expected to be an int fraction decimal or float xxx we should revisit whether using fractions to accumulate exact ratios is the right way to go the integer ratios for binary floats can have numerators or denominators with over 300 decimal digits the problem is more acute with decimal floats where the default decimal context supports a huge range of exponents from emin 999999 to emax 999999 when expanded with as_integer_ratio numbers like decimal 3 14e 5000 and decimal 3 14e 5000 have large numerators or denominators that will slow computation when the integer ratios are accumulated as fractions the size grows to cover the full range from the smallest magnitude to the largest for example fraction 3 14e 300 fraction 3 14e 300 has a 616 digit numerator likewise fraction decimal 3 14e 5000 fraction decimal 3 14e 5000 has 10 003 digit numerator this doesn t seem to have been problem in practice but it is a potential pitfall float nan or inf x may be an integral abc convert value to given numeric type t this covers the cases where t is fraction or where value is a nan or inf decimal or float fixme what do we do if this overflows iterate over values failing if any are less than zero rank order a dataset the lowest value has rank 1 ties are averaged so that equal values receive the same rank data 31 56 31 25 75 18 _rank data 3 5 5 0 3 5 2 0 6 0 1 0 the operation is idempotent _rank 3 5 5 0 3 5 2 0 6 0 1 0 3 5 5 0 3 5 2 0 6 0 1 0 it is possible to rank the data in reverse order so that the highest value has rank 1 also a key function can extract the field to be ranked goals eagles 45 bears 48 lions 44 _rank goals key itemgetter 1 reverse true 2 0 1 0 3 0 ranks are conventionally numbered starting from one however setting start to zero allows the ranks to be used as array indices prize gold silver bronze certificate scores 8 1 7 3 9 4 8 3 prize int i for i in _rank scores start 0 reverse true bronze certificate gold silver if this function becomes public at some point more thought needs to be given to the signature a list of ints is plausible when ties is min or max when ties is average either list float or list fraction is plausible default handling of ties matches scipy stats mstats spearmanr square root of n m rounded to the nearest integer using round to odd reference https www lri fr melquion doc 05 imacs17_1 expose pdf for 53 bit precision floats the bit width used in _float_sqrt_of_frac is 109 square root of n m as a float correctly rounded see principle and proof sketch at https bugs python org msg407078 convert to float square root of n m as a decimal correctly rounded premise for decimal computing n m sqrt can be off by 1 ulp from the correctly rounded result method check the result moving up or down a step if needed test n m root plus 2 2 test n m root minus 2 2 measures of central tendency averages return the sample arithmetic mean of data mean 1 2 3 4 4 2 8 from fractions import fraction as f mean f 3 7 f 1 21 f 5 3 f 1 3 fraction 13 21 from decimal import decimal as d mean d 0 5 d 0 75 d 0 625 d 0 375 decimal 0 5625 if data is empty statisticserror will be raised convert data to floats and compute the arithmetic mean this runs faster than the mean function and it always returns a float if the input dataset is empty it raises a statisticserror fmean 3 5 4 0 5 25 4 25 handle iterators that do not define __len__ convert data to floats and compute the geometric mean raises a statisticserror if the input dataset is empty if it contains a zero or if it contains a negative value no special efforts are made to achieve exact results however this may change in the future round geometric_mean 54 24 36 9 36 0 return the harmonic mean of data the harmonic mean is the reciprocal of the arithmetic mean of the reciprocals of the data it can be used for averaging ratios or rates for example speeds suppose a car travels 40 km hr for 5 km and then speeds up to 60 km hr for another 5 km what is the average speed harmonic_mean 40 60 48 0 suppose a car travels 40 km hr for 5 km and when traffic clears speeds up to 60 km hr for the remaining 30 km of the journey what is the average speed harmonic_mean 40 60 weights 5 30 56 0 if data is empty or any element is less than zero harmonic_mean will raise statisticserror fixme investigate ways to calculate medians without sorting quickselect return the median middle value of numeric data when the number of data points is odd return the middle data point when the number of data points is even the median is interpolated by taking the average of the two middle values median 1 3 5 3 median 1 3 5 7 4 0 return the low median of numeric data when the number of data points is odd the middle value is returned when it is even the smaller of the two middle values is returned median_low 1 3 5 3 median_low 1 3 5 7 3 return the high median of data when the number of data points is odd the middle value is returned when it is even the larger of the two middle values is returned median_high 1 3 5 3 median_high 1 3 5 7 5 estimates the median for numeric data binned around the midpoints of consecutive fixed width intervals the data can be any iterable of numeric data with each value being exactly the midpoint of a bin at least one value must be present the interval is width of each bin for example demographic information may have been summarized into consecutive ten year age groups with each group being represented by the 5 year midpoints of the intervals demographics counter 25 172 20 to 30 years old 35 484 30 to 40 years old 45 387 40 to 50 years old 55 22 50 to 60 years old 65 6 60 to 70 years old the 50th percentile median is the 536th person out of the 1071 member cohort that person is in the 30 to 40 year old age group the regular median function would assume that everyone in the tricenarian age group was exactly 35 years old a more tenable assumption is that the 484 members of that age group are evenly distributed between 30 and 40 for that we use median_grouped data list demographics elements median data 35 round median_grouped data interval 10 1 37 5 the caller is responsible for making sure the data points are separated by exact multiples of interval this is essential for getting a correct result the function does not check this precondition inputs may be any numeric type that can be coerced to a float during the interpolation step find the value at the midpoint remember this corresponds to the midpoint of the class interval using o log n bisection find where all the x values occur in the data all x will lie within data i j coerce to floats raising a typeerror if not possible interpolate the median using the formula found at https www cuemath com data median of grouped data lower limit of the median interval cumulative frequency of the preceding interval number of elements in the median internal return the most common data point from discrete or nominal data mode assumes discrete data and returns a single value this is the standard treatment of the mode as commonly taught in schools mode 1 1 2 3 3 3 3 4 3 this also works with nominal non numeric data mode red blue blue red green red red red if there are multiple modes with same frequency return the first one encountered mode red red green blue blue red if data is empty mode raises statisticserror return a list of the most frequently occurring values will return more than one result if there are multiple modes or an empty list if data is empty multimode aabbbbbbbbcc b multimode aabbbbccddddeeffffgg b d f multimode notes on methods for computing quantiles there is no one perfect way to compute quantiles here we offer two methods that serve common needs most other packages surveyed offered at least one or both of these two making them standard in the sense of widely adopted and reproducible they are also easy to explain easy to compute manually and have straight forward interpretations that aren t surprising the default method is known as r6 percentile exc or expected value of rank order statistics the alternative method is known as r7 percentile inc or mode of rank order statistics for sample data where there is a positive probability for values beyond the range of the data the r6 exclusive method is a reasonable choice consider a random sample of nine values from a population with a uniform distribution from 0 0 to 1 0 the distribution of the third ranked sample point is described by betavariate alpha 3 beta 7 which has mode 0 250 median 0 286 and mean 0 300 only the latter which corresponds with r6 gives the desired cut point with 30 of the population falling below that value making it comparable to a result from an inv_cdf function the r6 exclusive method is also idempotent for describing population data where the end points are known to be included in the data the r7 inclusive method is a reasonable choice instead of the mean it uses the mode of the beta distribution for the interior points per hyndman fan one nice property is that the vertices of q7 p divide the range into n 1 intervals and exactly 100p of the intervals lie to the left of q7 p and 100 1 p of the intervals lie to the right of q7 p if needed other methods could be added however for now the position is that fewer options make for easier choices and that external packages can be used for anything more advanced divide data into n continuous intervals with equal probability returns a list of n 1 cut points separating the intervals set n to 4 for quartiles the default set n to 10 for deciles set n to 100 for percentiles which gives the 99 cuts points that separate data in to 100 equal sized groups the data can be any iterable containing sample the cut points are linearly interpolated between data points if method is set to inclusive data is treated as population data the minimum value is treated as the 0th percentile and the maximum value is treated as the 100th percentile rescale i to m n clamp to 1 ld 1 exact integer math measures of spread see http mathworld wolfram com variance html http mathworld wolfram com samplevariance html return the sample variance of data data should be an iterable of real valued numbers with at least two values the optional argument xbar if given should be the mean of the data if it is missing or none the mean is automatically calculated use this function when your data is a sample from a population to calculate the variance from the entire population see pvariance examples data 2 75 1 75 1 25 0 25 0 5 1 25 3 5 variance data 1 3720238095238095 if you have already calculated the mean of your data you can pass it as the optional second argument xbar to avoid recalculating it m mean data variance data m 1 3720238095238095 this function does not check that xbar is actually the mean of data giving arbitrary values for xbar may lead to invalid or impossible results decimals and fractions are supported from decimal import decimal as d variance d 27 5 d 30 25 d 30 25 d 34 5 d 41 75 decimal 31 01875 from fractions import fraction as f variance f 1 6 f 1 2 f 5 3 fraction 67 108 return the population variance of data data should be a sequence or iterable of real valued numbers with at least one value the optional argument mu if given should be the mean of the data if it is missing or none the mean is automatically calculated use this function to calculate the variance from the entire population to estimate the variance from a sample the variance function is usually a better choice examples data 0 0 0 25 0 25 1 25 1 5 1 75 2 75 3 25 pvariance data 1 25 if you have already calculated the mean of the data you can pass it as the optional second argument to avoid recalculating it mu mean data pvariance data mu 1 25 decimals and fractions are supported from decimal import decimal as d pvariance d 27 5 d 30 25 d 30 25 d 34 5 d 41 75 decimal 24 815 from fractions import fraction as f pvariance f 1 4 f 5 4 f 1 2 fraction 13 72 return the square root of the sample variance see variance for arguments and other details stdev 1 5 2 5 2 5 2 75 3 25 4 75 1 0810874155219827 return the square root of the population variance see pvariance for arguments and other details pstdev 1 5 2 5 2 5 2 75 3 25 4 75 0 986893273527251 in one pass compute the mean and sample standard deviation as floats handle nans and infs gracefully finite inputs overflowed so scale down and recompute sqrt 1 sys float_info max non zero inputs underflowed so scale up and recompute scale 1 sqrt sys float_info min sys float_info epsilon improve accuracy with a differential correction https www wolframalpha com input i maclaurin series sqrt 28h 2 2b x 29 at x 3d0 statistics for relations between two inputs see https en wikipedia org wiki covariance https en wikipedia org wiki pearson_correlation_coefficient https en wikipedia org wiki simple_linear_regression covariance return the sample covariance of two inputs x and y covariance is a measure of the joint variability of two inputs x 1 2 3 4 5 6 7 8 9 y 1 2 3 1 2 3 1 2 3 covariance x y 0 75 z 9 8 7 6 5 4 3 2 1 covariance x z 7 5 covariance z x 7 5 pearson s correlation coefficient return the pearson s correlation coefficient for two inputs pearson s correlation coefficient r takes values between 1 and 1 it measures the strength and direction of a linear relationship x 1 2 3 4 5 6 7 8 9 y 9 8 7 6 5 4 3 2 1 correlation x x 1 0 correlation x y 1 0 if method is ranked computes spearman s rank correlation coefficient for two inputs the data is replaced by ranks ties are averaged so that equal values receive the same rank the resulting coefficient measures the strength of a monotonic relationship spearman s rank correlation coefficient is appropriate for ordinal data or for continuous data that doesn t meet the linear proportion requirement for pearson s correlation coefficient center rankings around zero slope and intercept for simple linear regression return the slope and intercept of simple linear regression parameters estimated using ordinary least squares simple linear regression describes relationship between an independent variable x and a dependent variable y in terms of a linear function y slope x intercept noise where slope and intercept are the regression parameters that are estimated and noise represents the variability of the data that was not explained by the linear regression it is equal to the difference between predicted and actual values of the dependent variable the parameters are returned as a named tuple x 1 2 3 4 5 noise normaldist samples 5 seed 42 y 3 x i 2 noise i for i in range 5 linear_regression x y doctest ellipsis linearregression slope 3 17495 intercept 1 00925 if proportional is true the independent variable x and the dependent variable y are assumed to be directly proportional the data is fit to a line passing through the origin since the intercept will always be 0 0 the underlying linear function simplifies to y slope x noise y 3 x i noise i for i in range 5 linear_regression x y proportional true doctest ellipsis linearregression slope 2 90475 intercept 0 0 list because used three times below generator because only used once below add zero to coerce result to a float equivalent to covariance x y variance x normal distribution there is no closed form solution to the inverse cdf for the normal distribution so we use a rational approximation instead wichura m j 1988 algorithm as241 the percentage points of the normal distribution applied statistics blackwell publishing 37 3 477 484 doi 10 2307 2347330 jstor 2347330 hash sum 55 88319_28806_14901_4439 hash sum 49 33206_50330_16102_89036 hash sum 47 52583_31754_92896_71629 if available use c implementation https en wikipedia org wiki normal_distribution https en wikipedia org wiki variance properties inverse cumulative distribution function x p x x p finds the value of the random variable such that the probability of the variable being less than or equal to that value equals the given probability this function is also called the percent point function or quantile function divide into n continuous intervals with equal probability returns a list of n 1 cut points separating the intervals set n to 4 for quartiles the default set n to 10 for deciles set n to 100 for percentiles which gives the 99 cuts points that separate the normal distribution in to 100 equal sized groups compute the overlapping coefficient ovl between two normal distributions measures the agreement between two normal probability distributions returns a value between 0 0 and 1 0 giving the overlapping area in the two underlying probability density functions n1 normaldist 2 4 1 6 n2 normaldist 3 2 2 0 n1 overlap n2 0 8035050657330205 see the overlapping coefficient as a measure of agreement between probability distributions and point estimation of the overlap of two normal densities henry f inman and edwin l bradley jr http dx doi org 10 1080 03610928908830127 sort to assure commutativity compute the standard score x mean stdev describes x in terms of the number of standard deviations above or below the mean of the normal distribution https www statisticshowto com probability and statistics z score return the mode of the normal distribution the mode is the value x where which the probability density function pdf takes its maximum value add a constant or another normaldist instance if other is a constant translate mu by the constant leaving sigma unchanged if other is a normaldist add both the means and the variances mathematically this works only if the two distributions are independent or if they are jointly normally distributed subtract a constant or another normaldist instance if other is a constant translate by the constant mu leaving sigma unchanged if other is a normaldist subtract the means and add the variances mathematically this works only if the two distributions are independent or if they are jointly normally distributed multiply both mu and sigma by a constant used for rescaling perhaps to change measurement units sigma is scaled with the absolute value of the constant divide both mu and sigma by a constant used for rescaling perhaps to change measurement units sigma is scaled with the absolute value of the constant
__all__ = [
    'NormalDist',
    'StatisticsError',
    'correlation',
    'covariance',
    'fmean',
    'geometric_mean',
    'harmonic_mean',
    'linear_regression',
    'mean',
    'median',
    'median_grouped',
    'median_high',
    'median_low',
    'mode',
    'multimode',
    'pstdev',
    'pvariance',
    'quantiles',
    'stdev',
    'variance',
]

import math
import numbers
import random
import sys

from fractions import Fraction
from decimal import Decimal
from itertools import count, groupby, repeat
from bisect import bisect_left, bisect_right
from math import hypot, sqrt, fabs, exp, erf, tau, log, fsum, sumprod
from math import isfinite, isinf
from functools import reduce
from operator import itemgetter
from collections import Counter, namedtuple, defaultdict

_SQRT2 = sqrt(2.0)



class StatisticsError(ValueError):
    pass




def _sum(data):
    
    count = 0
    types = set()
    types_add = types.add
    partials = {}
    partials_get = partials.get
    for typ, values in groupby(data, type):
        types_add(typ)
        for n, d in map(_exact_ratio, values):
            count += 1
            partials[d] = partials_get(d, 0) + n
    if None in partials:
        
        
        total = partials[None]
        assert not _isfinite(total)
    else:
        
        total = sum(Fraction(n, d) for d, n in partials.items())
    T = reduce(_coerce, types, int)  
    return (T, total, count)


def _ss(data, c=None):
    
    if c is not None:
        T, ssd, count = _sum((d := x - c) * d for x in data)
        return (T, ssd, c, count)
    count = 0
    types = set()
    types_add = types.add
    sx_partials = defaultdict(int)
    sxx_partials = defaultdict(int)
    for typ, values in groupby(data, type):
        types_add(typ)
        for n, d in map(_exact_ratio, values):
            count += 1
            sx_partials[d] += n
            sxx_partials[d] += n * n
    if not count:
        ssd = c = Fraction(0)
    elif None in sx_partials:
        
        
        ssd = c = sx_partials[None]
        assert not _isfinite(ssd)
    else:
        sx = sum(Fraction(n, d) for d, n in sx_partials.items())
        sxx = sum(Fraction(n, d*d) for d, n in sxx_partials.items())
        
        
        ssd = (count * sxx - sx * sx) / count
        c = sx / count
    T = reduce(_coerce, types, int)  
    return (T, ssd, c, count)


def _isfinite(x):
    try:
        return x.is_finite()  
    except AttributeError:
        return math.isfinite(x)  


def _coerce(T, S):
    
    
    assert T is not bool, "initial type T is bool"
    
    
    
    if T is S:  return T
    
    if S is int or S is bool:  return T
    if T is int:  return S
    
    if issubclass(S, T):  return S
    if issubclass(T, S):  return T
    
    if issubclass(T, int):  return S
    if issubclass(S, int):  return T
    
    if issubclass(T, Fraction) and issubclass(S, float):
        return S
    if issubclass(T, float) and issubclass(S, Fraction):
        return T
    
    msg = "don't know how to coerce %s and %s"
    raise TypeError(msg % (T.__name__, S.__name__))


def _exact_ratio(x):
    

    
    

    
    
    
    
    
    
    

    
    
    
    
    
    

    
    

    try:
        return x.as_integer_ratio()
    except AttributeError:
        pass
    except (OverflowError, ValueError):
        
        assert not _isfinite(x)
        return (x, None)
    try:
        
        return (x.numerator, x.denominator)
    except AttributeError:
        msg = f"can't convert type '{type(x).__name__}' to numerator/denominator"
        raise TypeError(msg)


def _convert(value, T):
    
    if type(value) is T:
        
        
        return value
    if issubclass(T, int) and value.denominator != 1:
        T = float
    try:
        
        return T(value)
    except TypeError:
        if issubclass(T, Decimal):
            return T(value.numerator) / T(value.denominator)
        else:
            raise


def _fail_neg(values, errmsg='negative value'):
    
    for x in values:
        if x < 0:
            raise StatisticsError(errmsg)
        yield x


def _rank(data, /, *, key=None, reverse=False, ties='average', start=1) -> list[float]:
    
    
    
    
    

    
    if ties != 'average':
        raise ValueError(f'Unknown tie resolution method: {ties!r}')
    if key is not None:
        data = map(key, data)
    val_pos = sorted(zip(data, count()), reverse=reverse)
    i = start - 1
    result = [0] * len(val_pos)
    for _, g in groupby(val_pos, key=itemgetter(0)):
        group = list(g)
        size = len(group)
        rank = i + (size + 1) / 2
        for value, orig_pos in group:
            result[orig_pos] = rank
        i += size
    return result


def _integer_sqrt_of_frac_rto(n: int, m: int) -> int:
    
    
    a = math.isqrt(n // m)
    return a | (a*a*m != n)




_sqrt_bit_width: int = 2 * sys.float_info.mant_dig + 3


def _float_sqrt_of_frac(n: int, m: int) -> float:
    
    
    q = (n.bit_length() - m.bit_length() - _sqrt_bit_width) // 2
    if q >= 0:
        numerator = _integer_sqrt_of_frac_rto(n, m << 2 * q) << q
        denominator = 1
    else:
        numerator = _integer_sqrt_of_frac_rto(n << -2 * q, m)
        denominator = 1 << -q
    return numerator / denominator   


def _decimal_sqrt_of_frac(n: int, m: int) -> Decimal:
    
    
    
    
    if n <= 0:
        if not n:
            return Decimal('0.0')
        n, m = -n, -m

    root = (Decimal(n) / Decimal(m)).sqrt()
    nr, dr = root.as_integer_ratio()

    plus = root.next_plus()
    np, dp = plus.as_integer_ratio()
    
    if 4 * n * (dr*dp)**2 > m * (dr*np + dp*nr)**2:
        return plus

    minus = root.next_minus()
    nm, dm = minus.as_integer_ratio()
    
    if 4 * n * (dr*dm)**2 < m * (dr*nm + dm*nr)**2:
        return minus

    return root




def mean(data):
    
    T, total, n = _sum(data)
    if n < 1:
        raise StatisticsError('mean requires at least one data point')
    return _convert(total / n, T)


def fmean(data, weights=None):
    
    if weights is None:
        try:
            n = len(data)
        except TypeError:
            
            n = 0
            def count(iterable):
                nonlocal n
                for n, x in enumerate(iterable, start=1):
                    yield x
            data = count(data)
        total = fsum(data)
        if not n:
            raise StatisticsError('fmean requires at least one data point')
        return total / n
    if not isinstance(weights, (list, tuple)):
        weights = list(weights)
    try:
        num = sumprod(data, weights)
    except ValueError:
        raise StatisticsError('data and weights must be the same length')
    den = fsum(weights)
    if not den:
        raise StatisticsError('sum of weights must be non-zero')
    return num / den


def geometric_mean(data):
    
    try:
        return exp(fmean(map(log, data)))
    except ValueError:
        raise StatisticsError('geometric mean requires a non-empty dataset '
                              'containing positive numbers') from None


def harmonic_mean(data, weights=None):
    
    if iter(data) is data:
        data = list(data)
    errmsg = 'harmonic mean does not support negative values'
    n = len(data)
    if n < 1:
        raise StatisticsError('harmonic_mean requires at least one data point')
    elif n == 1 and weights is None:
        x = data[0]
        if isinstance(x, (numbers.Real, Decimal)):
            if x < 0:
                raise StatisticsError(errmsg)
            return x
        else:
            raise TypeError('unsupported type')
    if weights is None:
        weights = repeat(1, n)
        sum_weights = n
    else:
        if iter(weights) is weights:
            weights = list(weights)
        if len(weights) != n:
            raise StatisticsError('Number of weights does not match data size')
        _, sum_weights, _ = _sum(w for w in _fail_neg(weights, errmsg))
    try:
        data = _fail_neg(data, errmsg)
        T, total, count = _sum(w / x if w else 0 for w, x in zip(weights, data))
    except ZeroDivisionError:
        return 0
    if total <= 0:
        raise StatisticsError('Weighted sum must be positive')
    return _convert(sum_weights / total, T)


def median(data):
    
    data = sorted(data)
    n = len(data)
    if n == 0:
        raise StatisticsError("no median for empty data")
    if n % 2 == 1:
        return data[n // 2]
    else:
        i = n // 2
        return (data[i - 1] + data[i]) / 2


def median_low(data):
    
    data = sorted(data)
    n = len(data)
    if n == 0:
        raise StatisticsError("no median for empty data")
    if n % 2 == 1:
        return data[n // 2]
    else:
        return data[n // 2 - 1]


def median_high(data):
    
    data = sorted(data)
    n = len(data)
    if n == 0:
        raise StatisticsError("no median for empty data")
    return data[n // 2]


def median_grouped(data, interval=1.0):
    
    data = sorted(data)
    n = len(data)
    if not n:
        raise StatisticsError("no median for empty data")

    
    
    x = data[n // 2]

    
    
    i = bisect_left(data, x)
    j = bisect_right(data, x, lo=i)

    
    try:
        interval = float(interval)
        x = float(x)
    except ValueError:
        raise TypeError(f'Value cannot be converted to a float')

    
    
    L = x - interval / 2.0    
    cf = i                    
    f = j - i                 
    return L + interval * (n / 2 - cf) / f


def mode(data):
    
    pairs = Counter(iter(data)).most_common(1)
    try:
        return pairs[0][0]
    except IndexError:
        raise StatisticsError('no mode for empty data') from None


def multimode(data):
    
    counts = Counter(iter(data))
    if not counts:
        return []
    maxcount = max(counts.values())
    return [value for value, count in counts.items() if count == maxcount]







































def quantiles(data, *, n=4, method='exclusive'):
    
    if n < 1:
        raise StatisticsError('n must be at least 1')
    data = sorted(data)
    ld = len(data)
    if ld < 2:
        if ld == 1:
            return data * (n - 1)
        raise StatisticsError('must have at least one data point')
    if method == 'inclusive':
        m = ld - 1
        result = []
        for i in range(1, n):
            j, delta = divmod(i * m, n)
            interpolated = (data[j] * (n - delta) + data[j + 1] * delta) / n
            result.append(interpolated)
        return result
    if method == 'exclusive':
        m = ld + 1
        result = []
        for i in range(1, n):
            j = i * m // n                               
            j = 1 if j < 1 else ld-1 if j > ld-1 else j  
            delta = i*m - j*n                            
            interpolated = (data[j - 1] * (n - delta) + data[j] * delta) / n
            result.append(interpolated)
        return result
    raise ValueError(f'Unknown method: {method!r}')








def variance(data, xbar=None):
    
    T, ss, c, n = _ss(data, xbar)
    if n < 2:
        raise StatisticsError('variance requires at least two data points')
    return _convert(ss / (n - 1), T)


def pvariance(data, mu=None):
    
    T, ss, c, n = _ss(data, mu)
    if n < 1:
        raise StatisticsError('pvariance requires at least one data point')
    return _convert(ss / n, T)


def stdev(data, xbar=None):
    
    T, ss, c, n = _ss(data, xbar)
    if n < 2:
        raise StatisticsError('stdev requires at least two data points')
    mss = ss / (n - 1)
    if issubclass(T, Decimal):
        return _decimal_sqrt_of_frac(mss.numerator, mss.denominator)
    return _float_sqrt_of_frac(mss.numerator, mss.denominator)


def pstdev(data, mu=None):
    
    T, ss, c, n = _ss(data, mu)
    if n < 1:
        raise StatisticsError('pstdev requires at least one data point')
    mss = ss / n
    if issubclass(T, Decimal):
        return _decimal_sqrt_of_frac(mss.numerator, mss.denominator)
    return _float_sqrt_of_frac(mss.numerator, mss.denominator)


def _mean_stdev(data):
    
    T, ss, xbar, n = _ss(data)
    if n < 2:
        raise StatisticsError('stdev requires at least two data points')
    mss = ss / (n - 1)
    try:
        return float(xbar), _float_sqrt_of_frac(mss.numerator, mss.denominator)
    except AttributeError:
        
        return float(xbar), float(xbar) / float(ss)

def _sqrtprod(x: float, y: float) -> float:
    "Return sqrt(x * y) computed with improved accuracy and without overflow/underflow."
    h = sqrt(x * y)
    if not isfinite(h):
        if isinf(h) and not isinf(x) and not isinf(y):
            
            scale = 2.0 ** -512  
            return _sqrtprod(scale * x, scale * y) / scale
        return h
    if not h:
        if x and y:
            
            
            scale = 2.0 ** 537
            return _sqrtprod(scale * x, scale * y) / scale
        return h
    
    
    d = sumprod((x, h), (y, -h))
    return h + d / (2.0 * h)









def covariance(x, y, /):
    
    n = len(x)
    if len(y) != n:
        raise StatisticsError('covariance requires that both inputs have same number of data points')
    if n < 2:
        raise StatisticsError('covariance requires at least two data points')
    xbar = fsum(x) / n
    ybar = fsum(y) / n
    sxy = sumprod((xi - xbar for xi in x), (yi - ybar for yi in y))
    return sxy / (n - 1)


def correlation(x, y, /, *, method='linear'):
    
    n = len(x)
    if len(y) != n:
        raise StatisticsError('correlation requires that both inputs have same number of data points')
    if n < 2:
        raise StatisticsError('correlation requires at least two data points')
    if method not in {'linear', 'ranked'}:
        raise ValueError(f'Unknown method: {method!r}')
    if method == 'ranked':
        start = (n - 1) / -2            
        x = _rank(x, start=start)
        y = _rank(y, start=start)
    else:
        xbar = fsum(x) / n
        ybar = fsum(y) / n
        x = [xi - xbar for xi in x]
        y = [yi - ybar for yi in y]
    sxy = sumprod(x, y)
    sxx = sumprod(x, x)
    syy = sumprod(y, y)
    try:
        return sxy / _sqrtprod(sxx, syy)
    except ZeroDivisionError:
        raise StatisticsError('at least one of the inputs is constant')


LinearRegression = namedtuple('LinearRegression', ('slope', 'intercept'))


def linear_regression(x, y, /, *, proportional=False):
    
    n = len(x)
    if len(y) != n:
        raise StatisticsError('linear regression requires that both inputs have same number of data points')
    if n < 2:
        raise StatisticsError('linear regression requires at least two data points')
    if not proportional:
        xbar = fsum(x) / n
        ybar = fsum(y) / n
        x = [xi - xbar for xi in x]  
        y = (yi - ybar for yi in y)  
    sxy = sumprod(x, y) + 0.0        
    sxx = sumprod(x, x)
    try:
        slope = sxy / sxx   
    except ZeroDivisionError:
        raise StatisticsError('x is constant')
    intercept = 0.0 if proportional else ybar - slope * xbar
    return LinearRegression(slope=slope, intercept=intercept)





def _normal_dist_inv_cdf(p, mu, sigma):
    
    
    
    
    
    q = p - 0.5
    if fabs(q) <= 0.425:
        r = 0.180625 - q * q
        
        num = (((((((2.50908_09287_30122_6727e+3 * r +
                     3.34305_75583_58812_8105e+4) * r +
                     6.72657_70927_00870_0853e+4) * r +
                     4.59219_53931_54987_1457e+4) * r +
                     1.37316_93765_50946_1125e+4) * r +
                     1.97159_09503_06551_4427e+3) * r +
                     1.33141_66789_17843_7745e+2) * r +
                     3.38713_28727_96366_6080e+0) * q
        den = (((((((5.22649_52788_52854_5610e+3 * r +
                     2.87290_85735_72194_2674e+4) * r +
                     3.93078_95800_09271_0610e+4) * r +
                     2.12137_94301_58659_5867e+4) * r +
                     5.39419_60214_24751_1077e+3) * r +
                     6.87187_00749_20579_0830e+2) * r +
                     4.23133_30701_60091_1252e+1) * r +
                     1.0)
        x = num / den
        return mu + (x * sigma)
    r = p if q <= 0.0 else 1.0 - p
    r = sqrt(-log(r))
    if r <= 5.0:
        r = r - 1.6
        
        num = (((((((7.74545_01427_83414_07640e-4 * r +
                     2.27238_44989_26918_45833e-2) * r +
                     2.41780_72517_74506_11770e-1) * r +
                     1.27045_82524_52368_38258e+0) * r +
                     3.64784_83247_63204_60504e+0) * r +
                     5.76949_72214_60691_40550e+0) * r +
                     4.63033_78461_56545_29590e+0) * r +
                     1.42343_71107_49683_57734e+0)
        den = (((((((1.05075_00716_44416_84324e-9 * r +
                     5.47593_80849_95344_94600e-4) * r +
                     1.51986_66563_61645_71966e-2) * r +
                     1.48103_97642_74800_74590e-1) * r +
                     6.89767_33498_51000_04550e-1) * r +
                     1.67638_48301_83803_84940e+0) * r +
                     2.05319_16266_37758_82187e+0) * r +
                     1.0)
    else:
        r = r - 5.0
        
        num = (((((((2.01033_43992_92288_13265e-7 * r +
                     2.71155_55687_43487_57815e-5) * r +
                     1.24266_09473_88078_43860e-3) * r +
                     2.65321_89526_57612_30930e-2) * r +
                     2.96560_57182_85048_91230e-1) * r +
                     1.78482_65399_17291_33580e+0) * r +
                     5.46378_49111_64114_36990e+0) * r +
                     6.65790_46435_01103_77720e+0)
        den = (((((((2.04426_31033_89939_78564e-15 * r +
                     1.42151_17583_16445_88870e-7) * r +
                     1.84631_83175_10054_68180e-5) * r +
                     7.86869_13114_56132_59100e-4) * r +
                     1.48753_61290_85061_48525e-2) * r +
                     1.36929_88092_27358_05310e-1) * r +
                     5.99832_20655_58879_37690e-1) * r +
                     1.0)
    x = num / den
    if q < 0.0:
        x = -x
    return mu + (x * sigma)



try:
    from _statistics import _normal_dist_inv_cdf
except ImportError:
    pass


class NormalDist:
    "Normal distribution of a random variable"
    
    

    __slots__ = {
        '_mu': 'Arithmetic mean of a normal distribution',
        '_sigma': 'Standard deviation of a normal distribution',
    }

    def __init__(self, mu=0.0, sigma=1.0):
        "NormalDist where mu is the mean and sigma is the standard deviation."
        if sigma < 0.0:
            raise StatisticsError('sigma must be non-negative')
        self._mu = float(mu)
        self._sigma = float(sigma)

    @classmethod
    def from_samples(cls, data):
        "Make a normal distribution instance from sample data."
        return cls(*_mean_stdev(data))

    def samples(self, n, *, seed=None):
        "Generate *n* samples for a given mean and standard deviation."
        rnd = random.random if seed is None else random.Random(seed).random
        inv_cdf = _normal_dist_inv_cdf
        mu = self._mu
        sigma = self._sigma
        return [inv_cdf(rnd(), mu, sigma) for _ in repeat(None, n)]

    def pdf(self, x):
        "Probability density function.  P(x <= X < x+dx) / dx"
        variance = self._sigma * self._sigma
        if not variance:
            raise StatisticsError('pdf() not defined when sigma is zero')
        diff = x - self._mu
        return exp(diff * diff / (-2.0 * variance)) / sqrt(tau * variance)

    def cdf(self, x):
        "Cumulative distribution function.  P(X <= x)"
        if not self._sigma:
            raise StatisticsError('cdf() not defined when sigma is zero')
        return 0.5 * (1.0 + erf((x - self._mu) / (self._sigma * _SQRT2)))

    def inv_cdf(self, p):
        
        if p <= 0.0 or p >= 1.0:
            raise StatisticsError('p must be in the range 0.0 < p < 1.0')
        return _normal_dist_inv_cdf(p, self._mu, self._sigma)

    def quantiles(self, n=4):
        
        return [self.inv_cdf(i / n) for i in range(1, n)]

    def overlap(self, other):
        
        
        
        
        
        if not isinstance(other, NormalDist):
            raise TypeError('Expected another NormalDist instance')
        X, Y = self, other
        if (Y._sigma, Y._mu) < (X._sigma, X._mu):  
            X, Y = Y, X
        X_var, Y_var = X.variance, Y.variance
        if not X_var or not Y_var:
            raise StatisticsError('overlap() not defined when sigma is zero')
        dv = Y_var - X_var
        dm = fabs(Y._mu - X._mu)
        if not dv:
            return 1.0 - erf(dm / (2.0 * X._sigma * _SQRT2))
        a = X._mu * Y_var - Y._mu * X_var
        b = X._sigma * Y._sigma * sqrt(dm * dm + dv * log(Y_var / X_var))
        x1 = (a + b) / dv
        x2 = (a - b) / dv
        return 1.0 - (fabs(Y.cdf(x1) - X.cdf(x1)) + fabs(Y.cdf(x2) - X.cdf(x2)))

    def zscore(self, x):
        
        
        if not self._sigma:
            raise StatisticsError('zscore() not defined when sigma is zero')
        return (x - self._mu) / self._sigma

    @property
    def mean(self):
        "Arithmetic mean of the normal distribution."
        return self._mu

    @property
    def median(self):
        "Return the median of the normal distribution"
        return self._mu

    @property
    def mode(self):
        
        return self._mu

    @property
    def stdev(self):
        "Standard deviation of the normal distribution."
        return self._sigma

    @property
    def variance(self):
        "Square of the standard deviation."
        return self._sigma * self._sigma

    def __add__(x1, x2):
        
        if isinstance(x2, NormalDist):
            return NormalDist(x1._mu + x2._mu, hypot(x1._sigma, x2._sigma))
        return NormalDist(x1._mu + x2, x1._sigma)

    def __sub__(x1, x2):
        
        if isinstance(x2, NormalDist):
            return NormalDist(x1._mu - x2._mu, hypot(x1._sigma, x2._sigma))
        return NormalDist(x1._mu - x2, x1._sigma)

    def __mul__(x1, x2):
        
        return NormalDist(x1._mu * x2, x1._sigma * fabs(x2))

    def __truediv__(x1, x2):
        
        return NormalDist(x1._mu / x2, x1._sigma / fabs(x2))

    def __pos__(x1):
        "Return a copy of the instance."
        return NormalDist(x1._mu, x1._sigma)

    def __neg__(x1):
        "Negates mu while keeping sigma the same."
        return NormalDist(-x1._mu, x1._sigma)

    __radd__ = __add__

    def __rsub__(x1, x2):
        "Subtract a NormalDist from a constant or another NormalDist."
        return -(x1 - x2)

    __rmul__ = __mul__

    def __eq__(x1, x2):
        "Two NormalDist objects are equal if their mu and sigma are both equal."
        if not isinstance(x2, NormalDist):
            return NotImplemented
        return x1._mu == x2._mu and x1._sigma == x2._sigma

    def __hash__(self):
        "NormalDist objects hash equal if their mu and sigma are both equal."
        return hash((self._mu, self._sigma))

    def __repr__(self):
        return f'{type(self).__name__}(mu={self._mu!r}, sigma={self._sigma!r})'

    def __getstate__(self):
        return self._mu, self._sigma

    def __setstate__(self, state):
        self._mu, self._sigma = state