Index: modules/math/ChangeLog ================================================================== --- modules/math/ChangeLog +++ modules/math/ChangeLog @@ -1,5 +1,21 @@ +2015-04-29 Arjen Markus + * statistics.tcl: Add test-Duckworth + * pdf_stat.tcl: Add empirical-distribution + * statistics.test: Add tests for test-Duckworth and empirical-distribution + * statistics.man: Describe test-Duckworth and empirical-distribution + +2015-04-28 Arjen Markus + * statistics.tcl: Bump version to 1.0 - Aku found the cause of earlier problems + +2015-04-26 Arjen Markus + * statistics.tcl: Bump version to 0.9.3 + Implemented an alternative to histogram (ticket 1502400fff) + Revised test-normal to use "significance" (ticket 2812473fff) + * statistics.man: Describe histogram-alt, changes to test-normal (and t-test-mean, again "confidence") + * pdf_stat.tcl: Correct the returned value for pdf-beta - if x is 0 or 1. + 2014-09-27 Arjen Markus * statistics.tcl: Bump version to 0.9.2 * statistics.test: Add tests for all pdf-* and cdf-* procedures, crude tests for random-* procedures * pdf_stat.tcl: Fix a typo (cdf-uniform) and fix inadvertent integer divisions should arguments be integer * special.tcl: Adding Christian's implementation of the inverse normal distribution function (invnorm) Index: modules/math/TODO ================================================================== --- modules/math/TODO +++ modules/math/TODO @@ -1,6 +1,15 @@ This file records outstanding actions for the math module + +dd. 26 april 2015, Arjen Markus +Add: +- additional linear algebra procedures by Federico Ferri +- lognormal income library by Eric Benedict +- empirical distribution +- tukey-duckworth test + + dd. 18 january 2014, Arjen Markus test cases for kernel-density: One test case is troublesome - uniform kernel, checking the total density Index: modules/math/calculus.tcl ================================================================== --- modules/math/calculus.tcl +++ modules/math/calculus.tcl @@ -10,11 +10,11 @@ # # RCS: @(#) $Id: calculus.tcl,v 1.15 2008/10/08 03:30:48 andreas_kupries Exp $ package require Tcl 8.4 package require math::interpolate -package provide math::calculus 0.8 +package provide math::calculus 0.8.1 # math::calculus -- # Namespace for the commands namespace eval ::math::calculus { @@ -1486,10 +1486,18 @@ 0.2044329400752989e+00 0.2094821410847278e+00} set qk15_wg { 0.1294849661688697e+00 0.2797053914892767e+00 0.3818300505051189e+00 0.4179591836734694e+00} } + +if {[package vsatisfies [package present Tcl] 8.5]} { + proc ::math::calculus::Min {a b} { expr {min ($a, $b)} } + proc ::math::calculus::Max {a b} { expr {max ($a, $b)} } +} else { + proc ::math::calculus::Min {a b} { if {$a < $b} { return $a } else { return $b }} + proc ::math::calculus::Max {a b} { if {$a > $b} { return $a } else { return $b }} +} proc ::math::calculus::qk15_basic {xstart xend func} { variable qk15_wg variable qk15_wgk variable qk15_xgk @@ -1555,14 +1563,14 @@ set result [expr {$resk*$hlgth}] set resabs [expr {$resabs*$dhlgth}] set resasc [expr {$resasc*$dhlgth}] set abserr [expr {abs(($resk-$resg)*$hlgth)}] if { $resasc != 0.0e+00 && $abserr != 0.0e+00 } { - set abserr [expr {$resasc*min(0.1e+01,(0.2e+3*$abserr/$resasc)**1.5e+00)}] + set abserr [expr {$resasc*[Min 0.1e+01 [expr {pow((0.2e+3*$abserr/$resasc),1.5e+00)}]]}] } if { $resabs > $uflow/(0.5e+02*$epmach) } { - set abserr [expr {max(($epmach*0.5e+02)*$resabs,$abserr)}] + set abserr [Max [expr {($epmach*0.5e+02)*$resabs}] $abserr] } return [list $result $abserr $resabs $resasc] } @@ -1623,15 +1631,15 @@ set resasc 0.0 for {set i 0} {$i < $n} {incr i} { set xb [expr {$xstart + $dx * $i}] set xe [expr {$xstart + $dx * ($i+1)}] - lassign [qk15_basic $xb $xe $func] dresult dabserr dresabs dresasc + foreach {dresult dabserr dresabs dresasc} [qk15_basic $xb $xe $func] break set result [expr {$result + $dresult}] set abserr [expr {$abserr + $dabserr}] set resabs [expr {$resabs + $dresabs}] set resasc [expr {$resasc + $dresasc}] } } return [list $result $abserr $resabs $resasc] } Index: modules/math/linalg.tcl ================================================================== --- modules/math/linalg.tcl +++ modules/math/linalg.tcl @@ -1740,11 +1740,11 @@ set U {} incr m -1 foreach row $A { lappend U [lrange $row 0 $m] } - puts $U + #puts $U } return [list $U $S $V] } # eigenvectorsSVD -- Index: modules/math/mvlinreg.tcl ================================================================== --- modules/math/mvlinreg.tcl +++ modules/math/mvlinreg.tcl @@ -3,11 +3,10 @@ # Copyright 2007 Eric Kemp-Benedict # Released under the BSD license under any terms # that allow it to be compatible with tcllib package require math::linearalgebra 1.0 -package require math::statistics 0.4 # ::math::statistics -- # This file adds: # mvlinreg = Multivariate Linear Regression # Index: modules/math/pdf_stat.tcl ================================================================== --- modules/math/pdf_stat.tcl +++ modules/math/pdf_stat.tcl @@ -12,22 +12,26 @@ # ::math::statistics -- # Namespace holding the procedures and variables # namespace eval ::math::statistics { - namespace export pdf-normal pdf-uniform \ + namespace export pdf-normal pdf-uniform pdf-lognormal \ pdf-exponential \ - cdf-normal cdf-uniform \ + cdf-normal cdf-uniform cdf-lognormal \ cdf-exponential \ cdf-students-t \ - random-normal random-uniform \ + random-normal random-uniform random-lognormal \ random-exponential \ histogram-uniform \ pdf-gamma pdf-poisson pdf-chisquare pdf-students-t pdf-beta \ + pdf-weibull pdf-gumbel pdf-pareto pdf-cauchy \ cdf-gamma cdf-poisson cdf-chisquare cdf-beta \ + cdf-weibull cdf-gumbel cdf-pareto cdf-cauchy \ random-gamma random-poisson random-chisquare random-students-t random-beta \ - incompleteGamma incompleteBeta + random-weibull random-gumbel random-pareto random-cauchy \ + incompleteGamma incompleteBeta \ + estimate-pareto empirical-distribution variable cdf_normal_prob {} variable cdf_normal_x {} variable cdf_toms322_cached {} variable initialised_cdf 0 @@ -56,10 +60,40 @@ } set xn [expr {($x-$mean)/$stdev}] set prob [expr {exp(-$xn*$xn/2.0)/$stdev/$factorNormalPdf}] + return $prob +} + +# pdf-lognormal -- +# Return the probabilities belonging to a log-normal distribution +# +# Arguments: +# mean Mean of the distribution +# stdev Standard deviation +# x Value for which the probability must be determined +# +# Result: +# Probability of value x under the given distribution +# +proc ::math::statistics::pdf-lognormal { mean stdev x } { + variable NEGSTDEV + variable factorNormalPdf + + if { $stdev <= 0.0 || $mean <= 0.0 } { + return -code error -errorcode ARG \ + -errorinfo "Standard deviation and mean must be positive" \ + "Standard deviation and mean must be positive" + } + + set sigma [expr {sqrt(log(1.0 + $stdev /double($mean*$mean)))}] + set mu [expr {log($mean) - 0.5 * $sigma * $sigma}] + + set xn [expr {(log($x)-$mu)/$sigma}] + set prob [expr {exp(-$xn*$xn/2.0)/$sigma/$factorNormalPdf}] + return $prob } # pdf-uniform -- @@ -139,10 +173,45 @@ set xn [expr {($x-double($mean))/$stdev}] set prob1 [Cdf-toms322 1 5000 [expr {$xn*$xn}]] if { $xn > 0.0 } { set prob [expr {0.5+0.5*$prob1}] + } else { + set prob [expr {0.5-0.5*$prob1}] + } + + return $prob +} + + +# cdf-lognormal -- +# Return the cumulative probability belonging to a log-normal distribution +# +# Arguments: +# mean Mean of the distribution +# stdev Standard deviation +# x Value for which the probability must be determined +# +# Result: +# Probability of value x under the given distribution +# +proc ::math::statistics::cdf-lognormal { mean stdev x } { + variable NEGSTDEV + + if { $stdev <= 0.0 || $mean <= 0.0 } { + return -code error -errorcode ARG \ + -errorinfo "Standard deviation and mean must be positive" \ + "Standard deviation and mean must be positive" + } + + set sigma [expr {sqrt(log(1.0 + $stdev /double($mean*$mean)))}] + set mu [expr {log($mean) - 0.5 * $sigma * $sigma}] + + set xn [expr {(log($x)-$mu)/$sigma}] + set prob1 [Cdf-toms322 1 5000 [expr {$xn*$xn}]] + if { $xn > 0.0 } { + set prob [expr {0.5+0.5*$prob1}] } else { set prob [expr {0.5-0.5*$prob1}] } return $prob @@ -493,10 +562,63 @@ return $result } + +# random-lognormal -- +# Return a list of random numbers satisfying a log-normal +# distribution +# +# Arguments: +# mean Mean of the distribution +# stdev Standard deviation of the distribution +# number Number of values to generate +# +# Result: +# List of random numbers +# +# Note: +# This version uses the Box-Muller transformation, +# a quick and robust method for generating normally- +# distributed numbers. +# +proc ::math::statistics::random-lognormal { mean stdev number } { + variable twopi + + if { $stdev <= 0.0 || $mean <= 0.0 } { + return -code error -errorcode ARG \ + -errorinfo "Standard deviation and mean must be positive" \ + "Standard deviation and mean must be positive" + } + + set sigma [expr {sqrt(log(1.0 + $stdev /double($mean*$mean)))}] + set mu [expr {log($mean) - 0.5 * $sigma * $sigma}] + +# set result {} +# for { set i 0 } {$i < $number } { incr i } { +# lappend result [Inverse-cdf-normal $mean $stdev [expr {rand()}]] +# } + + puts "Random-lognormal: $mu -- $sigma" + + set result {} + + for { set i 0 } {$i < $number } { incr i 2 } { + set angle [expr {$twopi * rand()}] + set rad [expr {sqrt(-2.0*log(rand()))}] + set xrand [expr {$rad * cos($angle)}] + set yrand [expr {$rad * sin($angle)}] + lappend result [expr {exp($mu + $sigma * $xrand)}] + if { $i < $number-1 } { + lappend result [expr {exp($mu + $sigma * $yrand)}] + } + } + + return $result +} + # Cdf-toms322 -- # Calculate the cumulative density function for several distributions # according to TOMS322 # # Arguments: @@ -971,14 +1093,14 @@ } # # Corner cases ... need to check these! # if {$x == 0.0} { - return 0.0 + return [expr {$a > 1.0? 0.0 : Inf}] } if {$x == 1.0} { - return 0.0 + return [expr {$b > 1.0? 0.0 : Inf}] } set aplusb [expr {$a + $b}] set term1 [expr {[::math::ln_Gamma $aplusb]- [::math::ln_Gamma $a] - [::math::ln_Gamma $b]}] set term2 [expr {($a - 1.0) * log($x) + ($b - 1.0) * log(1.0 - $x)}] @@ -1109,10 +1231,111 @@ } return $retval } + +# pdf-weibull -- +# Return the probabilities belonging to a Weibull distribution +# +# Arguments: +# scale Scale parameter of the Weibull distribution +# shape Shape parameter of the Weibull distribution +# x Value of variate +# +# Result: +# Probability density of the given value of x to occur +# +# Note: +# "-$x ** $shape" is evaluated as "(-$x)**$shape", hence use a division +# +proc ::math::statistics::pdf-weibull { scale shape x } { + variable OUTOFRANGE + + if { $x < 0 } { + return 0.0 + } + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {$x / double($scale)}] + return [expr {$shape/double($scale) * pow($x,($shape-1.0)) / exp(pow($x,$shape))}] +} + + +# pdf-gumbel -- +# Return the probabilities belonging to a Gumbel distribution +# +# Arguments: +# location Location parameter of the Gumbel distribution +# scale Scale parameter of the Gumbel distribution +# x Value of variate +# +# Result: +# Probability density of the given value of x to occur +# +proc ::math::statistics::pdf-gumbel { location scale x } { + variable OUTOFRANGE + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {($x - $location) / double($scale)}] + return [expr {exp(-$x - exp(-$x)) / $scale}] +} + + +# pdf-pareto -- +# Return the probabilities belonging to a Pareto distribution +# +# Arguments: +# scale Scale parameter of the Pareto distribution +# shape Shape parameter of the Pareto distribution +# x Value of variate +# +# Result: +# Probability density of the given value of x to occur +# +proc ::math::statistics::pdf-pareto { scale shape x } { + variable OUTOFRANGE + + if { $x <= $scale } { + return 0.0 + } + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {$x / double($scale)}] + return [expr {$shape / double($scale) / pow($x,($shape + 1.0))}] +} + + +# pdf-cauchy -- +# Return the probabilities belonging to a Cauchy distribution +# +# Arguments: +# location Location parameter of the Cauchy distribution +# scale Scale parameter of the Cauchy distribution +# x Value of variate +# +# Result: +# Probability density of the given value of x to occur +# +# Note: +# The Cauchy distribution does not have finite higher-order moments +# +proc ::math::statistics::pdf-cauchy { location scale x } { + variable OUTOFRANGE + variable pi + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {($x - $location) / double($scale)}] + return [expr {1.0 / $pi / $scale / (1.0 +$x*$x)}] +} + # cdf-gamma -- # Return the cumulative probabilities belonging to a gamma distribution # # Arguments: @@ -1192,18 +1415,113 @@ # proc ::math::statistics::cdf-beta { a b x } { incompleteBeta $a $b $x } + +# cdf-weibull -- +# Return the cumulative probabilities belonging to a Weibull distribution +# +# Arguments: +# scale Scale parameter of the Weibull distribution +# shape Shape parameter of the Weibull distribution +# x Value of variate +# +# Result: +# Cumulative probability of the given value of x to occur +# +proc ::math::statistics::cdf-weibull { scale shape x } { + variable OUTOFRANGE + + if { $x <= 0 } { + return 0.0 + } + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {$x / double($scale)}] + return [expr {1.0 - 1.0 / exp(pow($x,$shape))}] +} + + +# cdf-gumbel -- +# Return the cumulative probabilities belonging to a Gumbel distribution +# +# Arguments: +# location Location parameter of the Gumbel distribution +# scale Scale parameter of the Gumbel distribution +# x Value of variate +# +# Result: +# Cumulative probability of the given value of x to occur +# +proc ::math::statistics::cdf-gumbel { location scale x } { + variable OUTOFRANGE + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {($x - $location) / double($scale)}] + return [expr {exp( -exp(-$x) )}] +} + + +# cdf-pareto -- +# Return the cumulative probabilities belonging to a Pareto distribution +# +# Arguments: +# scale Scale parameter of the Pareto distribution +# shape Shape parameter of the Pareto distribution +# x Value of variate +# +# Result: +# Cumulative probability density of the given value of x to occur +# +proc ::math::statistics::cdf-pareto { scale shape x } { + variable OUTOFRANGE + + if { $x <= $scale } { + return 0.0 + } + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {$x / double($scale)}] + return [expr {1.0 - 1.0 / pow($x,$shape)}] +} + + +# cdf-cauchy -- +# Return the cumulative probabilities belonging to a Cauchy distribution +# +# Arguments: +# location Scale parameter of the Cauchy distribution +# scale Shape parameter of the Cauchy distribution +# x Value of variate +# +# Result: +# Cumulative probability density of the given value of x to occur +# +proc ::math::statistics::cdf-cauchy { location scale x } { + variable OUTOFRANGE + variable pi + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set x [expr {($x - $location) / double($scale)}] + return [expr {0.5 + atan($x) / $pi}] +} + # random-gamma -- # Generate a list of gamma-distributed deviates # # Arguments: # alpha Shape parameter # beta Rate parameter -# x Value of variate +# number Number of values to return # # Result: # List of random values # # Note: @@ -1459,10 +1777,206 @@ lappend retval $k } return $retval } + +# random-weibull -- +# Generate a list of Weibull distributed deviates +# +# Arguments: +# scale Scale parameter of the Weibull distribution +# shape Shape parameter of the Weibull distribution +# number Number of values to return +# +# Result: +# List of random values +# +proc ::math::statistics::random-weibull { scale shape number } { + variable OUTOFRANGE + + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set rshape [expr {1.0/$shape}] + set retval {} + for {set i 0} {$i < $number} {incr i} { + lappend retval [expr {$scale * pow( (-log(rand())),$rshape)}] + } + return $retval +} + + +# random-gumbel -- +# Generate a list of Weibull distributed deviates +# +# Arguments: +# location Location parameter of the Gumbel distribution +# scale Scale parameter of the Gumbel distribution +# number Number of values to return +# +# Result: +# List of random values +# +proc ::math::statistics::random-gumbel { location scale number } { + variable OUTOFRANGE + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set retval {} + for {set i 0} {$i < $number} {incr i} { + lappend retval [expr {$location - $scale * log(-log(rand()))}] + } + return $retval +} + + +# random-pareto -- +# Generate a list of Pareto distributed deviates +# +# Arguments: +# scale Scale parameter of the Pareto distribution +# shape Shape parameter of the Pareto distribution +# number Number of values to return +# +# Result: +# List of random values +# +proc ::math::statistics::random-pareto { scale shape number } { + variable OUTOFRANGE + + if { $scale <= 0.0 || $shape <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set rshape [expr {1.0/$shape}] + set retval {} + for {set i 0} {$i < $number} {incr i} { + lappend retval [expr {$scale / pow(rand(),$rshape)}] + } + return $retval +} + + +# random-cauchy -- +# Generate a list of Cauchy distributed deviates +# +# Arguments: +# location Location parameter of the Cauchy distribution +# scale Shape parameter of the Cauchy distribution +# number Number of values to return +# +# Result: +# List of random values +# +proc ::math::statistics::random-cauchy { location scale number } { + variable OUTOFRANGE + variable pi + + if { $scale <= 0.0 } { + return -code error -errorcode ARG -errorinfo $OUTOFRANGE $OUTOFRANGE + } + set retval {} + for {set i 0} {$i < $number} {incr i} { + lappend retval [expr {$location + $scale * tan( $pi * (rand() - 0.5))}] + } + return $retval +} + + +# estimate-pareto -- +# Estimate the parameters of a Pareto distribution +# +# Arguments: +# values Values that are supposed to be distributed according to Pareto +# +# Result: +# Estimates of the scale and shape parameters as well as the standard error +# for the shape parameter. +# +proc ::math::statistics::estimate-pareto { values } { + variable OUTOFRANGE + variable TOOFEWDATA + + set nvalues {} + set negative 0 + + foreach v $values { + if { $v != {} } { + lappend nvalues $v + if { $v <= 0.0 } { + set negative 1 + } + } + } + if { [llength $nvalues] == 0 } { + return -code error -errorcode ARG -errorinfo $TOOFEWDATA $TOOFEWDATA + } + if { $negative } { + return -code error -errorcode ARG -errorinfo "One or more negative or zero values" $OUTOFRANGE + } + + # + # Scale parameter + # + set scale [min $nvalues] + + # + # Shape parameter + # + set n [llength $nvalues] + set sum 0.0 + foreach v $nvalues { + set sum [expr {$sum + log($v) - log($scale)}] + } + set shape [expr {$n / $sum}] + + return [list $scale $shape [expr {$shape/sqrt($n)}]] +} + + +# empirical-distribution -- +# Determine the empirical distribution +# +# Arguments: +# values Values that are to be examined +# +# Result: +# List of sorted values and their empirical probability +# +# Note: +# The value of "a" is adopted from the corresponding Wikipedia page, +# which in turn adopted it from the R "stats" package (qqnorm function) +# +proc ::math::statistics::empirical-distribution { values } { + variable TOOFEWDATA + + set n [llength $values] + + if { $n < 5 } { + return -code error -errorcode ARG -errorinfo $TOOFEWDATA $TOOFEWDATA + } + + set a 0.375 + if { $n > 10 } { + set a 0.5 + } + + set distribution {} + set idx 1 + foreach x [lsort -real -increasing $values] { + if { $x != {} } { + set p [expr {($idx - $a) / ($n + 1 - 2.0 * $a)}] + + lappend distribution $x $p + incr idx + } + } + + return $distribution +} + # # Simple numerical tests # if { [info exists ::argv0] && ([file tail [info script]] == [file tail $::argv0]) } { Index: modules/math/pkgIndex.tcl ================================================================== --- modules/math/pkgIndex.tcl +++ modules/math/pkgIndex.tcl @@ -12,13 +12,13 @@ if {![package vsatisfies [package provide Tcl] 8.3]} {return} package ifneeded math::roman 1.0 [list source [file join $dir romannumerals.tcl]] if {![package vsatisfies [package provide Tcl] 8.4]} {return} # statistics depends on linearalgebra (for multi-variate linear regression). -package ifneeded math::statistics 0.9.2 [list source [file join $dir statistics.tcl]] +package ifneeded math::statistics 0.9.3 [list source [file join $dir statistics.tcl]] package ifneeded math::optimize 1.0.1 [list source [file join $dir optimize.tcl]] -package ifneeded math::calculus 0.8 [list source [file join $dir calculus.tcl]] +package ifneeded math::calculus 0.8.1 [list source [file join $dir calculus.tcl]] package ifneeded math::interpolate 1.1 [list source [file join $dir interpolate.tcl]] package ifneeded math::linearalgebra 1.1.5 [list source [file join $dir linalg.tcl]] package ifneeded math::bignum 3.1.1 [list source [file join $dir bignum.tcl]] package ifneeded math::bigfloat 1.2.2 [list source [file join $dir bigfloat.tcl]] package ifneeded math::machineparameters 0.1 [list source [file join $dir machineparameters.tcl]] Index: modules/math/statistics.man ================================================================== --- modules/math/statistics.man +++ modules/math/statistics.man @@ -157,10 +157,22 @@ intervals of the histogram. [arg_def list values] - List of data [arg_def list weights] - List of weights, one weight per value [list_end] [para] + +[call [cmd ::math::statistics::histogram-alt] [arg limits] [arg values] [opt weights]] +Alternative implementation of the histogram procedure: the open end of the intervals +is at the lower bound instead of the upper bound. + +[list_begin arguments] +[arg_def list limits] - List of upper limits (in ascending order) for the +intervals of the histogram. +[arg_def list values] - List of data +[arg_def list weights] - List of weights, one weight per value +[list_end] +[para] [call [cmd ::math::statistics::corr] [arg data1] [arg data2]] Determine the correlation coefficient between two sets of data. [list_begin arguments] @@ -177,38 +189,45 @@ [list_begin arguments] [arg_def list data] - List of raw data values (small sample) [arg_def float confidence] - Confidence level (0.95 or 0.99 for instance) [list_end] [para] + [call [cmd ::math::statistics::t-test-mean] [arg data] [arg est_mean] \ -[arg est_stdev] [arg confidence]] +[arg est_stdev] [arg alpha]] Test whether the mean value of a sample is in accordance with the -estimated normal distribution with a certain level of confidence. +estimated normal distribution with a certain probability. Returns 1 if the test succeeds or 0 if the mean is unlikely to fit the given distribution. [list_begin arguments] [arg_def list data] - List of raw data values (small sample) [arg_def float est_mean] - Estimated mean of the distribution [arg_def float est_stdev] - Estimated stdev of the distribution -[arg_def float confidence] - Confidence level (0.95 or 0.99 for instance) +[arg_def float alpha] - Probability level (0.95 or 0.99 for instance) [list_end] [para] -[call [cmd ::math::statistics::test-normal] [arg data] [arg confidence]] + +[call [cmd ::math::statistics::test-normal] [arg data] [arg significance]] Test whether the given data follow a normal distribution -with a certain level of confidence. +with a certain level of significance. Returns 1 if the data are normally distributed within the level of -confidence, returns 0 if not. The underlying test is the Lilliefors -test. +significance, returns 0 if not. The underlying test is the Lilliefors +test. Smaller values of the significance mean a stricter testing. [list_begin arguments] [arg_def list data] - List of raw data values -[arg_def float confidence] - Confidence level (one of 0.80, 0.90, 0.95 or 0.99) +[arg_def float significance] - Significance level (one of 0.01, 0.05, 0.10, 0.15 or 0.20). For compatibility +reasons the values "1-significance", 0.80, 0.85, 0.90, 0.95 or 0.99 are also accepted. [list_end] [para] +Compatibility issue: the original implementation and documentation used the term "confidence" and used a value +1-significance (see ticket 2812473fff). This has been corrected as of version 0.9.3. +[para] + [call [cmd ::math::statistics::lillieforsFit] [arg data]] Returns the goodness of fit to a normal distribution according to Lilliefors. The higher the number, the more likely the data are indeed normally distributed. The test requires at least [emph five] data @@ -217,10 +236,36 @@ [list_begin arguments] [arg_def list data] - List of raw data values [list_end] [para] + +[call [cmd ::math::statistics::test-Duckworth] [arg list1] [arg list2] [arg significance]] +Determine if two data sets have the same median according to the Tukey-Duckworth test. +The procedure returns 0 if the medians are unequal, 1 if they are equal, -1 if the test can not +be conducted (the smallest value must be in a different set than the greatest value). +# +# Arguments: +# list1 Values in the first data set +# list2 Values in the second data set +# significance Significance level (either 0.05, 0.01 or 0.001) +# +# Returns: + +Test whether the given data follow a normal distribution +with a certain level of significance. +Returns 1 if the data are normally distributed within the level of +significance, returns 0 if not. The underlying test is the Lilliefors +test. Smaller values of the significance mean a stricter testing. + +[list_begin arguments] +[arg_def list list1] - First list of data +[arg_def list list2] - Second list of data +[arg_def float significance] - Significance level (either 0.05, 0.01 or 0.001) +[list_end] +[para] + [call [cmd ::math::statistics::quantiles] [arg data] [arg confidence]] Return the quantiles for a given set of data [list_begin arguments] [para] [arg_def list data] - List of raw data values @@ -643,25 +688,35 @@ [section "STATISTICAL DISTRIBUTIONS"] In the literature a large number of probability distributions can be found. The statistics package supports: [list_begin itemized] [item] -The normal or Gaussian distribution +The normal or Gaussian distribution as well as the log-normal distribution [item] The uniform distribution - equal probability for all data within a given interval [item] The exponential distribution - useful as a model for certain extreme-value distributions. [item] The gamma distribution - based on the incomplete Gamma integral +[item] +The beta distribution [item] The chi-square distribution [item] The student's T distribution [item] The Poisson distribution +[item] +The Pareto distribution +[item] +The Gumbel distribution +[item] +The Weibull distribution +[item] +The Cauchy distribution [item] PM - binomial,F. [list_end] In principle for each distribution one has procedures for: @@ -684,10 +739,21 @@ [call [cmd ::math::statistics::pdf-normal] [arg mean] [arg stdev] [arg value]] Return the probability of a given value for a normal distribution with given mean and standard deviation. +[list_begin arguments] +[arg_def float mean] - Mean value of the distribution +[arg_def float stdev] - Standard deviation of the distribution +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::pdf-lognormal] [arg mean] [arg stdev] [arg value]] +Return the probability of a given value for a log-normal distribution with +given mean and standard deviation. + [list_begin arguments] [arg_def float mean] - Mean value of the distribution [arg_def float stdev] - Standard deviation of the distribution [arg_def float value] - Value for which the probability is required [list_end] @@ -752,27 +818,95 @@ [list_begin arguments] [arg_def float df] - Degrees of freedom [arg_def float value] - Value for which the probability is required [list_end] [para] + +[call [cmd ::math::statistics::pdf-gamma] [arg a] [arg b] [arg value]] +Return the probability of a given value for a Gamma +distribution with given shape and rate parameters + +[list_begin arguments] +[arg_def float a] - Shape parameter +[arg_def float b] - Rate parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] [call [cmd ::math::statistics::pdf-beta] [arg a] [arg b] [arg value]] Return the probability of a given value for a Beta distribution with given shape parameters [list_begin arguments] [arg_def float a] - First shape parameter -[arg_def float b] - First shape parameter +[arg_def float b] - Second shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::pdf-weibull] [arg scale] [arg shape] [arg value]] +Return the probability of a given value for a Weibull +distribution with given scale and shape parameters + +[list_begin arguments] +[arg_def float location] - Scale parameter +[arg_def float scale] - Shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::pdf-gumbel] [arg location] [arg scale] [arg value]] +Return the probability of a given value for a Gumbel +distribution with given location and shape parameters + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::pdf-pareto] [arg scale] [arg shape] [arg value]] +Return the probability of a given value for a Pareto +distribution with given scale and shape parameters + +[list_begin arguments] +[arg_def float scale] - Scale parameter +[arg_def float shape] - Shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::pdf-cauchy] [arg location] [arg scale] [arg value]] +Return the probability of a given value for a Cauchy +distribution with given location and shape parameters. Note that the Cauchy distribution +has no finite higher-order moments. + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Shape parameter [arg_def float value] - Value for which the probability is required [list_end] [para] [call [cmd ::math::statistics::cdf-normal] [arg mean] [arg stdev] [arg value]] Return the cumulative probability of a given value for a normal distribution with given mean and standard deviation, that is the probability for values up to the given one. +[list_begin arguments] +[arg_def float mean] - Mean value of the distribution +[arg_def float stdev] - Standard deviation of the distribution +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::cdf-lognormal] [arg mean] [arg stdev] [arg value]] +Return the cumulative probability of a given value for a log-normal +distribution with given mean and standard deviation, that is the +probability for values up to the given one. + [list_begin arguments] [arg_def float mean] - Mean value of the distribution [arg_def float stdev] - Standard deviation of the distribution [arg_def float value] - Value for which the probability is required [list_end] @@ -808,11 +942,11 @@ [list_end] [para] [call [cmd ::math::statistics::cdf-gamma] [arg alpha] [arg beta] [arg value]] Return the cumulative probability of a given value for a Gamma -distribution with given shape and rate parameters +distribution with given shape and rate parameters. [list_begin arguments] [arg_def float alpha] - Shape parameter [arg_def float beta] - Rate parameter [arg_def float value] - Value for which the cumulative probability is required @@ -819,11 +953,11 @@ [list_end] [para] [call [cmd ::math::statistics::cdf-poisson] [arg mu] [arg k]] Return the cumulative probability of a given number of occurrences in -the same interval (k) for a Poisson distribution with given mean (mu) +the same interval (k) for a Poisson distribution with given mean (mu). [list_begin arguments] [arg_def float mu] - Mean number of occurrences [arg_def int k] - Number of occurences [list_end] @@ -833,18 +967,81 @@ Return the cumulative probability of a given value for a Beta distribution with given shape parameters [list_begin arguments] [arg_def float a] - First shape parameter -[arg_def float b] - First shape parameter +[arg_def float b] - Second shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::cdf-weibull] [arg scale] [arg shape] [arg value]] +Return the cumulative probability of a given value for a Weibull +distribution with given scale and shape parameters. + +[list_begin arguments] +[arg_def float scale] - Scale parameter +[arg_def float shape] - Shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::cdf-gumbel] [arg location] [arg scale] [arg value]] +Return the cumulative probability of a given value for a Gumbel +distribution with given location and scale parameters. + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Scale parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::cdf-pareto] [arg scale] [arg shape] [arg value]] +Return the cumulative probability of a given value for a Pareto +distribution with given scale and shape parameters + +[list_begin arguments] +[arg_def float scale] - Scale parameter +[arg_def float shape] - Shape parameter +[arg_def float value] - Value for which the probability is required +[list_end] +[para] + +[call [cmd ::math::statistics::cdf-cauchy] [arg location] [arg scale] [arg value]] +Return the cumulative probability of a given value for a Cauchy +distribution with given location and scale parameters. + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Scale parameter [arg_def float value] - Value for which the probability is required [list_end] [para] + +[call [cmd ::math::statistics::empirical-distribution] [arg values]] +Return a list of values and their empirical probability. The values are sorted in increasing order. +(The implementation follows the description at the corresponding Wikipedia page) + +[list_begin arguments] +[arg_def list values] - List of data to be examined +[list_end] +[para] [call [cmd ::math::statistics::random-normal] [arg mean] [arg stdev] [arg number]] Return a list of "number" random values satisfying a normal distribution with given mean and standard deviation. +[list_begin arguments] +[arg_def float mean] - Mean value of the distribution +[arg_def float stdev] - Standard deviation of the distribution +[arg_def int number] - Number of values to be returned +[list_end] +[para] + +[call [cmd ::math::statistics::random-lognormal] [arg mean] [arg stdev] [arg number]] +Return a list of "number" random values satisfying a log-normal +distribution with given mean and standard deviation. [list_begin arguments] [arg_def float mean] - Mean value of the distribution [arg_def float stdev] - Standard deviation of the distribution [arg_def int number] - Number of values to be returned [list_end] @@ -870,11 +1067,11 @@ [list_end] [para] [call [cmd ::math::statistics::random-gamma] [arg alpha] [arg beta] [arg number]] Return a list of "number" random values satisfying -a Gamma distribution with given shape and rate parameters +a Gamma distribution with given shape and rate parameters. [list_begin arguments] [arg_def float alpha] - Shape parameter [arg_def float beta] - Rate parameter [arg_def int number] - Number of values to be returned @@ -881,48 +1078,92 @@ [list_end] [para] [call [cmd ::math::statistics::random-poisson] [arg mu] [arg number]] Return a list of "number" random values satisfying -a Poisson distribution with given mean +a Poisson distribution with given mean. [list_begin arguments] [arg_def float mu] - Mean of the distribution [arg_def int number] - Number of values to be returned [list_end] [para] [call [cmd ::math::statistics::random-chisquare] [arg df] [arg number]] Return a list of "number" random values satisfying -a chi square distribution with given degrees of freedom +a chi square distribution with given degrees of freedom. [list_begin arguments] [arg_def float df] - Degrees of freedom [arg_def int number] - Number of values to be returned [list_end] [para] [call [cmd ::math::statistics::random-student-t] [arg df] [arg number]] Return a list of "number" random values satisfying -a Student's t distribution with given degrees of freedom +a Student's t distribution with given degrees of freedom. [list_begin arguments] [arg_def float df] - Degrees of freedom [arg_def int number] - Number of values to be returned [list_end] [para] [call [cmd ::math::statistics::random-beta] [arg a] [arg b] [arg number]] Return a list of "number" random values satisfying -a Beta distribution with given shape parameters +a Beta distribution with given shape parameters. [list_begin arguments] [arg_def float a] - First shape parameter [arg_def float b] - Second shape parameter [arg_def int number] - Number of values to be returned [list_end] [para] + +[call [cmd ::math::statistics::random-weibull] [arg scale] [arg shape] [arg number]] +Return a list of "number" random values satisfying +a Weibull distribution with given scale and shape parameters. + +[list_begin arguments] +[arg_def float scale] - Scale parameter +[arg_def float shape] - Shape parameter +[arg_def int number] - Number of values to be returned +[list_end] +[para] + +[call [cmd ::math::statistics::random-gumbel] [arg location] [arg scale] [arg number]] +Return a list of "number" random values satisfying +a Gumbel distribution with given location and scale parameters. + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Scale parameter +[arg_def int number] - Number of values to be returned +[list_end] +[para] + +[call [cmd ::math::statistics::random-pareto] [arg scale] [arg shape] [arg number]] +Return a list of "number" random values satisfying +a Pareto distribution with given scale and shape parameters. + +[list_begin arguments] +[arg_def float scale] - Scale parameter +[arg_def float shape] - Shape parameter +[arg_def int number] - Number of values to be returned +[list_end] +[para] + +[call [cmd ::math::statistics::random-cauchy] [arg location] [arg scale] [arg number]] +Return a list of "number" random values satisfying +a Cauchy distribution with given location and scale parameters. + +[list_begin arguments] +[arg_def float location] - Location parameter +[arg_def float scale] - Scale parameter +[arg_def int number] - Number of values to be returned +[list_end] +[para] [call [cmd ::math::statistics::histogram-uniform] [arg xmin] [arg xmax] [arg limits] [arg number]] Return the expected histogram for a uniform distribution. [list_begin arguments] @@ -957,10 +1198,19 @@ [arg_def float b] - Second shape parameter [arg_def float x] - Value of x (limit of the integral) [arg_def float tol] - Required tolerance (default: 1.0e-9) [list_end] [para] + +[call [cmd ::math::statistics::estimate-pareto] [arg values]] +Estimate the parameters for the Pareto distribution that comes closest to the given values. +Returns the estimated scale and shape parameters, as well as the standard error for the shape parameter. + +[list_begin arguments] +[arg_def list values] - List of values, assumed to be distributed according to a Pareto distribution +[list_end] +[para] [list_end] TO DO: more function descriptions to be added [section "DATA MANIPULATION"] Index: modules/math/statistics.tcl ================================================================== --- modules/math/statistics.tcl +++ modules/math/statistics.tcl @@ -15,13 +15,14 @@ # version 0.6: added pdf and cdf procedures for various distributions # (provided by Eric Kemp-Benedict) # version 0.7: added Kruskal-Wallis test (by Torsten Berg) # version 0.8: added Wilcoxon test and Spearman rank correlation # version 0.9: added kernel density estimation +# version 0.9.3: added histogram-alt, corrected test-normal package require Tcl 8.4 -package provide math::statistics 0.9.2 +package provide math::statistics 1.0 package require math if {![llength [info commands ::lrepeat]]} { # Forward portability, emulate lrepeat proc ::lrepeat {n args} { @@ -44,17 +45,18 @@ namespace eval ::math::statistics { # # Safer: change to short procedures # namespace export mean min max number var stdev pvar pstdev basic-stats corr \ - histogram interval-mean-stdev t-test-mean quantiles \ + histogram histogram-alt interval-mean-stdev t-test-mean quantiles \ test-normal lillieforsFit \ autocorr crosscorr filter map samplescount median \ test-2x2 print-2x2 control-xbar test_xbar \ control-Rchart test-Rchart \ test-Kruskal-Wallis analyse-Kruskal-Wallis group-rank \ - test-Wilcoxon spearman-rank spearman-rank-extended + test-Wilcoxon spearman-rank spearman-rank-extended \ + test-Duckworth # # Error messages # variable NEGSTDEV {Zero or negative standard deviation} variable TOOFEWDATA {Too few or invalid data} @@ -227,10 +229,75 @@ set index 0 set found 0 foreach limit $limits { if { $value <= $limit } { + set found 1 + set buckets($index) [expr $buckets($index)+$weight] + break + } + incr index + } + + if { $found == 0 } { + set buckets($last) [expr $buckets($last)+$weight] + } + } + + set result {} + for { set index 0 } { $index <= $last } { incr index } { + lappend result $buckets($index) + } + + return $result +} + +# histogram-alt -- +# Return histogram information from a list of numbers - +# intervals are open-ended at the lower bound instead of at the upper bound +# +# Arguments: +# limits Upper limits for the buckets (in increasing order) +# values List of values to be examined +# weights List of weights, one per value (optional) +# +# Results: +# List of number of values in each bucket (length is one more than +# the number of limits) +# +# +proc ::math::statistics::histogram-alt { limits values {weights {}} } { + + if { [llength $limits] < 1 } { + return -code error -errorcode ARG -errorinfo {No limits given} {No limits given} + } + if { [llength $weights] > 0 && [llength $values] != [llength $weights] } { + return -code error -errorcode ARG -errorinfo {Number of weights be equal to number of values} {Weights and values differ in length} + } + + set limits [lsort -real -increasing $limits] + + for { set index 0 } { $index <= [llength $limits] } { incr index } { + set buckets($index) 0 + } + + set last [llength $limits] + + # Will do integer arithmetic if unset + if {$weights eq ""} { + set weights [lrepeat [llength $values] 1] + } + + foreach value $values weight $weights { + if { $value == {} } { + continue + } + + set index 0 + set found 0 + foreach limit $limits { + if { $value < $limit } { set found 1 set buckets($index) [expr $buckets($index)+$weight] break } incr index @@ -430,67 +497,76 @@ # test-normal -- # Test for normality (using method Lilliefors) # # Arguments: # data Values that need to be tested -# confidence ... +# significance Level at which the discrepancy from normality is tested # # Result: -# 1 if of the statistic D +# 1 if the Lilliefors statistic D is larger than the critical level +# +# Note: +# There was a mistake in the implementation before 0.9.3: confidence (wrong word) +# instead of significance. To keep compatibility with earlier versions, both +# significance and 1-significance are accepted. # -proc ::math::statistics::test-normal {data confidence} { +proc ::math::statistics::test-normal {data significance} { set D [lillieforsFit $data] + + if { $significance > 0.5 } { + set significance [expr {1.0-$significance}] ;# Convert the erroneous levels pre 0.9.3 + } set Dcrit -- - if { abs($confidence-0.80) < 0.0001 } { + if { abs($significance-0.20) < 0.0001 } { set Dcrit 0.741 } - if { abs($confidence-0.85) < 0.0001 } { + if { abs($significance-0.15) < 0.0001 } { set Dcrit 0.775 } - if { abs($confidence-0.90) < 0.0001 } { + if { abs($significance-0.10) < 0.0001 } { set Dcrit 0.819 } - if { abs($confidence-0.95) < 0.0001 } { + if { abs($significance-0.05) < 0.0001 } { set Dcrit 0.895 } - if { abs($confidence-0.99) < 0.0001 } { + if { abs($significance-0.01) < 0.0001 } { set Dcrit 1.035 } if { $Dcrit != "--" } { return [expr {$D > $Dcrit ? 1 : 0 }] } else { - return -code error "Confidence level must be one of: 0.80, 0.85, 0.90, 0.95 or 0.99" + return -code error "Significancce level must be one of: 0.20, 0.15, 0.10, 0.05 or 0.01" } } # t-test-mean -- # Test whether the mean value of a sample is in accordance with the -# estimated normal distribution with a certain level of confidence +# estimated normal distribution with a certain probability # (Student's t test) # # Arguments: # data List of raw data values (small sample) # est_mean Estimated mean of the distribution # est_stdev Estimated stdev of the distribution -# confidence Confidence level (0.95 or 0.99 for instance) +# alpha Probability level (0.95 or 0.99 for instance) # # Result: # 1 if the test is positive, 0 otherwise. If there are too few data, # returns an empty string # -proc ::math::statistics::t-test-mean { data est_mean est_stdev confidence } { +proc ::math::statistics::t-test-mean { data est_mean est_stdev alpha } { variable NEGSTDEV variable TOOFEWDATA if { $est_stdev <= 0.0 } { return -code error -errorcode ARG -errorinfo $NEGSTDEV $NEGSTDEV } set allstats [BasicStats all $data] - set conf2 [expr {(1.0+$confidence)/2.0}] + set alpha2 [expr {(1.0+$alpha)/2.0}] set sample_mean [lindex $allstats 0] set sample_number [lindex $allstats 3] if { $sample_number > 1 } { @@ -497,11 +573,11 @@ set tzero [expr {abs($sample_mean-$est_mean)/$est_stdev * \ sqrt($sample_number-1)}] set degrees [expr {$sample_number-1}] set prob [cdf-students-t $degrees $tzero] - return [expr {$prob<$conf2}] + return [expr {$prob<$alpha2}] } else { return -code error -errorcode DATA -errorinfo $TOOFEWDATA $TOOFEWDATA } } @@ -1311,10 +1387,102 @@ if { $range > $rupper } { lappend result $i } } return $result } + +# test-Duckworth -- +# Determine if two data sets have the same median according to the Tukey-Duckworth test +# +# Arguments: +# list1 Values in the first data set +# list2 Values in the second data set +# significance Significance level (either 0.05, 0.01 or 0.001) +# +# Returns: +# 0 if the medians are unequal, 1 if they are equal, -1 if the test can not +# be conducted (the smallest value must be in a different set than the greatest value) +# +proc ::math::statistics::test-Duckworth {list1 list2 significance} { + set sorted1 [lsort -real $list1] + set sorted2 [lsort -real -decreasing $list2] + + set lowest1 [lindex $sorted1 0] + set lowest2 [lindex $sorted2 end] + set greatest1 [lindex $sorted1 end] + set greatest2 [lindex $sorted2 0] + + if { $lowest1 <= $lowest2 && $greatest1 >= $greatest2 } { + return -1 + } + if { $lowest1 >= $lowest2 && $greatest1 <= $greatest2 } { + return -1 + } + + # + # Determine how many elements of set 1 are lower than the lowest of set 2 + # Ditto for the number of elements of set 2 greater than the greatest of set 1 + # (Or vice versa) + # + if { $lowest1 < $lowest2 } { + set lowest $lowest2 + set greatest $greatest1 + } else { + set lowest $lowest1 + set greatest $greatest2 + set sorted1 [lsort -real $list2] + set sorted2 [lsort -real -decreasing $list1] + #lassign [list $sorted1 $sorted2] sorted2 sorted1 + } + + set count1 0 + set count2 0 + foreach v1 $sorted1 { + if { $v1 >= $lowest } { + break + } + incr count1 + } + foreach v2 $sorted2 { + if { $v2 <= $greatest } { + break + } + incr count2 + } + + # + # Determine the statistic D, possibly with correction + # + set n1 [llength $list1] + set n2 [llength $list2] + + set correction 0 + if { 3 + 4*$n1/3 <= $n2 && $n2 <= 2*$n1 } { + set correction -1 + } + if { 3 + 4*$n2/3 <= $n1 && $n1 <= 2*$n2 } { + set correction -1 + } + + set D [expr {$count1 + $count2 + $correction}] + + switch -- [string trim $significance 0] { + ".05" { + return [expr {$D >= 7? 0 : 1}] + } + ".01" { + return [expr {$D >= 10? 0 : 1}] + } + ".001" { + return [expr {$D >= 13? 0 : 1}] + } + default { + return -code error "Significance level must be 0.05, 0.01 or 0.001" + } + } +} + # # Load the auxiliary scripts # source [file join [file dirname [info script]] pdf_stat.tcl] Index: modules/math/statistics.test ================================================================== --- modules/math/statistics.test +++ modules/math/statistics.test @@ -198,10 +198,19 @@ test "Histogram-1.5" "Histogram - linear data 2 with weights" -match glob -body { set values [::math::statistics::histogram {1.5 2.5} [concat $::data_linear 0.0 0.0] \ [concat [lrepeat [llength $::data_linear] 1] 0 0]] } -result {1 1 8} +# +# Alternative definition of the intervals (ticket 1502400fff) +# Note the difference in the expected bin sizes for the two +# +test "Histogram-2.1" "Histogram - alternative interval bounds" -match glob -body { + set values [concat [::math::statistics::histogram-alt {5.0 7.0} $::data_linear] \ + [::math::statistics::histogram {5.0 7.0} $::data_linear]] +} -result {4 2 4 5 2 3} + # # Quantiles # Bug #1272910: related to rounding 0.5 - use different levels instead # because another bug was fixed, return to the original # levels again @@ -438,21 +447,43 @@ ::math::statistics::lillieforsFit {72 112 111 107 119 92 126 80 81 84 115 118 128 128 123 116 125 126 122 126 127 86 142 132 87 123 133 106 103 118 114 94} } -result 0.82827415657 -test "Testnormal-1.0" "Test birthweight data for normality - 80%" -match tolerant -body { +test "Testnormal-1.0" "Test birthweight data for normality - 20% significance" -match exact -body { ::math::statistics::test-normal {72 112 111 107 119 92 126 80 81 84 115 118 128 128 123 116 125 126 122 126 127 86 - 142 132 87 123 133 106 103 118 114 94} 0.80 + 142 132 87 123 133 106 103 118 114 94} 0.20 } -result 1 -test "Testnormal-1.0" "Test birthweight data for normality - 95%" -match tolerant -body { +test "Testnormal-1.0" "Test birthweight data for normality - 5% significance" -match exact -body { ::math::statistics::test-normal {72 112 111 107 119 92 126 80 81 84 115 118 128 128 123 116 125 126 122 126 127 86 - 142 132 87 123 133 106 103 118 114 94} 0.95 + 142 132 87 123 133 106 103 118 114 94} 0.05 } -result 0 + +test "Test-Duckworth-1.0" "Test Tukey-Duckworth - 5% significance" -match exact -body { + set list1 {10 2 3 4 6} + set list2 {12 3 4 6} + + ::math::statistics::test-Duckworth $list1 $list2 0.05 +} -result 1 + +test "Test-Duckworth-1.1" "Test Tukey-Duckworth - symmetry" -match exact -body { + set list1 {1 2 3 4 5 6 7 8 9 10} + set list2 {6 7 8 9 10 11 12 13 14 15 16 17} + + set result [list [::math::statistics::test-Duckworth $list1 $list2 0.05] \ + [::math::statistics::test-Duckworth $list2 $list1 0.05]] +} -result {0 0} + +test "Test-Duckworth-1.2" "Test Tukey-Duckworth - applicability" -match exact -body { + set list1 {2 3 4 6 20} + set list2 {12 3 4 6} + + ::math::statistics::test-Duckworth $list1 $list2 0.05 +} -result -1 # # Testing multivariate linear regression # # Provide some data @@ -537,10 +568,38 @@ [::math::statistics::cdf-normal 2.0 2.0 4.0] \ [::math::statistics::cdf-normal -2.0 2.0 0.0] \ [::math::statistics::cdf-normal 2.0 2.0 3.0]] } -result {0.8413205502059895 0.8413205502059895 0.8413205502059895 0.8413205502059895 0.691451459572962} +test "lognormal-distribution-1.0" "Test pdf-lognormal" -match tolerant -body { + foreach {mu sigma mean stdev} {0.0 1.0 mean1 stdev1 2.0 2.0 mean2 stdev2 -2.0 2.0 mean3 stdev3} { + set m [expr {exp($mu + $sigma*$sigma/2.0)}] + set $mean $m + set $stdev [expr {(exp($sigma*$sigma) - 1.0) * $m*$m}] + } + + set x [list \ + [::math::statistics::pdf-lognormal $mean1 $stdev1 [expr {exp(1.0)}]] \ + [::math::statistics::pdf-lognormal $mean2 $stdev2 [expr {exp(4.0)}]] \ + [::math::statistics::pdf-lognormal $mean3 $stdev3 [expr {exp(0.0)}]] \ + [::math::statistics::pdf-lognormal $mean2 $stdev2 [expr {exp(3.0)}]]] +} -result {0.24197072451914337 0.12098536225957168 0.12098536225957168 0.17603266338214976} + +test "lognormal-distribution-1.1" "Test cdf-lognormal" -match tolerant -body { + foreach {mu sigma mean stdev} {0.0 1.0 mean1 stdev1 2.0 2.0 mean2 stdev2 -2.0 2.0 mean3 stdev3} { + set m [expr {exp($mu + $sigma*$sigma/2.0)}] + set $mean $m + set $stdev [expr {(exp($sigma*$sigma) - 1.0) * $m*$m}] + } + + set x [list \ + [::math::statistics::cdf-lognormal $mean1 $stdev1 [expr {exp(1.0)}]] \ + [::math::statistics::cdf-lognormal $mean2 $stdev2 [expr {exp(4.0)}]] \ + [::math::statistics::cdf-lognormal $mean3 $stdev3 [expr {exp(0.0)}]] \ + [::math::statistics::cdf-lognormal $mean2 $stdev2 [expr {exp(3.0)}]]] +} -result {0.8413205502059895 0.8413205502059895 0.8413205502059895 0.691451459572962} + test "gamma-distribution-1.0" "Test pdf-gamma" -match tolerant -body { set x [list \ [::math::statistics::pdf-gamma 1.5 2.7 3.0] \ [::math::statistics::pdf-gamma 7.5 0.2 30.0] \ [::math::statistics::pdf-gamma 15.0 1.2 2.0]] @@ -625,39 +684,182 @@ } -result {0.16220409275804 0.998630771123192 1.0 0.000125234318666948 0.0728881294218269 2.99872547567313e-23 3.07056696205524e-09 0.998641008671625 0.765865005703006 0.999999999996 0.000125237075575121 8.23161135486914e-20 0.464369443974288 0.5 1.0 0.8208} +# +# TODO: chose the tests with _integer_ arguments more carefully +# +test "gumbel-distribution-1.0" "Test pdf-gumbel" -match tolerant -body { + set x [list \ + [::math::statistics::pdf-gumbel 1.0 1.0 0.0] \ + [::math::statistics::pdf-gumbel 1.0 1.0 0.1] \ + [::math::statistics::pdf-gumbel 1.0 1.0 0.2] \ + [::math::statistics::pdf-gumbel 1.0 1.0 1.0] \ + [::math::statistics::pdf-gumbel 1.0 1.0 2.0] \ + [::math::statistics::pdf-gumbel 1.0 1.0 5.0] \ + [::math::statistics::pdf-gumbel 0.1 2.0 0.0] \ + [::math::statistics::pdf-gumbel 0.1 2.0 1.0] \ + [::math::statistics::pdf-gumbel 0.1 2.0 2.0] \ + [::math::statistics::pdf-gumbel 0.1 2.0 5.0] \ + [::math::statistics::pdf-gumbel 1 1 5 ] ] +} -result {0.179374 0.210219 0.240378 0.367879 0.254646 0.017983 0.183706 0.168507 0.131350 0.039580 0.017983} + +test "gumbel-distribution-1.1" "Test cdf-gumbel" -match tolerant -body { + set x [list \ + [::math::statistics::cdf-gumbel 1.0 1.0 0.0] \ + [::math::statistics::cdf-gumbel 1.0 1.0 0.2] \ + [::math::statistics::cdf-gumbel 1.0 1.0 1.0] \ + [::math::statistics::cdf-gumbel 1.0 1.0 2.0] \ + [::math::statistics::cdf-gumbel 0.1 2.0 0.0] \ + [::math::statistics::cdf-gumbel 0.1 2.0 1.0] \ + [::math::statistics::cdf-gumbel 0.1 2.0 2.0] \ + [::math::statistics::cdf-gumbel 1 1 2 ] ] +} -result {0.065988 0.108009 0.367879 0.692201 0.349493 0.528544 0.679266 0.692201} + +test "weibull-distribution-1.0" "Test pdf-weibull" -match tolerant -body { + set x [list \ + [::math::statistics::pdf-weibull 1.0 1.0 -1.0] \ + [::math::statistics::pdf-weibull 1.0 1.0 0.0] \ + [::math::statistics::pdf-weibull 1.0 1.0 0.1] \ + [::math::statistics::pdf-weibull 1.0 1.0 0.2] \ + [::math::statistics::pdf-weibull 1.0 1.0 1.0] \ + [::math::statistics::pdf-weibull 1.0 1.0 2.0] \ + [::math::statistics::pdf-weibull 1.0 1.0 5.0] \ + [::math::statistics::pdf-weibull 2.0 2.0 0.0] \ + [::math::statistics::pdf-weibull 2.0 2.0 1.0] \ + [::math::statistics::pdf-weibull 2.0 2.0 2.0] \ + [::math::statistics::pdf-weibull 2.0 2.0 5.0] ] +} -result {0 1.0 0.904837 0.818730 0.367879 0.135335 0.006738 0 0.389400 0.367879 0.004826} + +test "weibull-distribution-1.1" "Test cdf-weibull" -match tolerant -body { + set x [list \ + [::math::statistics::cdf-weibull 1.0 1.0 -1.0] \ + [::math::statistics::cdf-weibull 1.0 1.0 0.0] \ + [::math::statistics::cdf-weibull 1.0 1.0 0.2] \ + [::math::statistics::cdf-weibull 1.0 1.0 1.0] \ + [::math::statistics::cdf-weibull 1.0 1.0 2.0] \ + [::math::statistics::cdf-weibull 2.0 2.0 0.0] \ + [::math::statistics::cdf-weibull 2.0 2.0 1.0] \ + [::math::statistics::cdf-weibull 2.0 2.0 2.0] \ + [::math::statistics::cdf-weibull 2 2 2 ] ] +} -result {0 0 0.181269 0.632106 0.864665 0 0.221199 0.632121 0.632121} + +test "pareto-distribution-1.0" "Test pdf-pareto" -match tolerant -body { + set x [list \ + [::math::statistics::pdf-pareto 1.0 1.0 0.0] \ + [::math::statistics::pdf-pareto 1.0 1.0 1.1] \ + [::math::statistics::pdf-pareto 1.0 1.0 1.2] \ + [::math::statistics::pdf-pareto 1.0 1.0 2.0] \ + [::math::statistics::pdf-pareto 1.0 1.0 3.0] \ + [::math::statistics::pdf-pareto 1.0 1.0 5.0] \ + [::math::statistics::pdf-pareto 2.0 2.0 2.1] \ + [::math::statistics::pdf-pareto 2.0 2.0 3.0] \ + [::math::statistics::pdf-pareto 2.0 2.0 5.0] \ + [::math::statistics::pdf-pareto 2.0 2.0 10.0] ] +} -result {0 0.826446 0.694444 0.25 0.111111 0.04 0.863838 0.296296 0.064 0.008} + +test "pareto-distribution-1.1" "Test cdf-pareto" -match tolerant -body { + set x [list \ + [::math::statistics::cdf-pareto 1.0 1.0 0.0] \ + [::math::statistics::cdf-pareto 1.0 1.0 1.1] \ + [::math::statistics::cdf-pareto 1.0 1.0 1.2] \ + [::math::statistics::cdf-pareto 1.0 1.0 2.0] \ + [::math::statistics::cdf-pareto 1.0 1.0 3.0] \ + [::math::statistics::cdf-pareto 2.0 2.0 2.1] \ + [::math::statistics::cdf-pareto 2.0 2.0 3.0] \ + [::math::statistics::cdf-pareto 2.0 2.0 5.0] \ + [::math::statistics::cdf-pareto 2 2 3 ] ] +} -result {0 0.090909 0.1666667 0.5 0.666667 0.092971 0.555556 0.84 0.555556} + +test "cauchy-distribution-1.0" "Test pdf-cauchy" -match tolerant -body { + set x [list \ + [::math::statistics::pdf-cauchy 1.0 1.0 0.0] \ + [::math::statistics::pdf-cauchy 2.0 1.0 1.0] \ + [::math::statistics::pdf-cauchy 1.0 2.0 2.0] \ + [::math::statistics::pdf-cauchy 2.0 2.0 2.0] ] +} -result {0.1591555 0.1591555 0.1273240 0.1591550} + +test "cauchy-distribution-1.1" "Test cdf-cauchy" -match tolerant -body { + set x [list \ + [::math::statistics::cdf-cauchy 1.0 1.0 0.0] \ + [::math::statistics::cdf-cauchy 2.0 1.0 1.0] \ + [::math::statistics::cdf-cauchy 1.0 2.0 2.0] \ + [::math::statistics::cdf-cauchy 2.0 2.0 2.0] ] +} -result {0.25 0.25 0.6475836 0.5} + +test "empirical-distribution-1.0" "Test empirical-distribution" -match tolerant -body { + set x {10 4 3 2 5 6 7} + set distribution [::math::statistics::empirical-distribution $x] +} -result {2 0.086207 3 0.224138 4 0.36207 5 0.5 6 0.637910 7 0.775862 10 0.913793} + # # Crude tests for the random number generators # Mainly to verify that there are no obvious errors # -test "random-numbers-1.0" "Test random-uniform" -body { - set rnumbers [::math::statistics::random-uniform 0 10 100] - - set result 1 - foreach r $rnumbers { - if { $r < 0.0 || $r > 10.0 } { - result 0 - break - } - } - lappend result [llength $rnumbers] -} -result {1 100} - -test "random-numbers-1.1" "Test random-exponential" -body { - set rnumbers [::math::statistics::random-exponential 1 100] - - set result 1 - foreach r $rnumbers { - if { $r < 0.0 } { - result 0 - break - } - } - lappend result [llength $rnumbers] -} -result {1 100} +# To verify that the values are scaled properly, use a fixed seed +# +set ::rseed 1000000 + +test "random-numbers-1.0" "Test random-uniform" -body { + expr {srand($::rseed)} + + set rnumbers [::math::statistics::random-uniform 0 10 100] + + set inrange 1 + foreach r $rnumbers { + if { $r < 0.0 || $r > 10.0 } { + set inrange 0 + break + } + } + + expr {srand($::rseed)} + set scaled 1 + set rnumbers2 [::math::statistics::random-uniform 0 20 100] + foreach r1 $rnumbers r2 $rnumbers2 { + set scale [expr {$r2 / $r1}] + if { abs($scale - 2.0) > 0.00001 } { + set scaled 0 + } + } + expr {srand($::rseed)} + set shifted 1 + set rnumbers3 [::math::statistics::random-uniform 10 20 100] + foreach r1 $rnumbers r3 $rnumbers3 { + set shift [expr {$r3 - $r1}] + if { abs($shift - 10.0) > 0.00001 } { + set shifted 0 + } + } + + set result [list $inrange [llength $rnumbers] $scaled $shifted] +} -result {1 100 1 1} + +test "random-numbers-1.1" "Test random-exponential" -body { + expr {srand($::rseed)} + set rnumbers [::math::statistics::random-exponential 1 100] + + set inrange 1 + foreach r $rnumbers { + if { $r < 0.0 } { + set inrange 0 + break + } + } + + expr {srand($::rseed)} + set scaled 1 + set rnumbers2 [::math::statistics::random-exponential 2 100] + foreach r1 $rnumbers r2 $rnumbers2 { + set scale [expr {$r2 / $r1}] + if { abs($scale - 2.0) > 0.00001 } { + set scaled 0 + } + } + set result [list $inrange [llength $rnumbers] $scaled] +} -result {1 100 1} test "random-numbers-1.2" "Test random-normal" -body { set rnumbers [::math::statistics::random-normal 0 1 100] set result [llength $rnumbers] } -result 100 @@ -691,10 +893,48 @@ break } } lappend result [llength $rnumbers] } -result {1 100} + +test "random-numbers-1.8" "Test random-gumbel" -body { + set rnumbers [::math::statistics::random-gumbel 1.0 3.0 100] + set result [llength $rnumbers] +} -result 100 + +test "random-numbers-1.9" "Test random-weibull" -body { + set rnumbers [::math::statistics::random-weibull 1.0 3.0 100] + set result [llength $rnumbers] +} -result 100 + +test "random-numbers-1.10" "Test random-pareto" -body { + set rnumbers [::math::statistics::random-pareto 1.0 3.0 100] + set result 1 + foreach r $rnumbers { + if { $r < 1.0 } { + result 0 + break + } + } + lappend result [llength $rnumbers] +} -result {1 100} + +test "random-numbers-1.11" "Test random-lognormal" -body { + set rnumbers [::math::statistics::random-lognormal 1 1 100] + set result [llength $rnumbers] +} -result 100 + +test "random-numbers-1.11" "Test random-cauchy" -body { + set rnumbers [::math::statistics::random-cauchy 0 1 100] + set result [llength $rnumbers] +} -result 100 + +test "random-numbers-2.1" "Test estimate-pareto" -match tolerant -body { + expr {srand($::rseed)} + set rnumbers [::math::statistics::random-pareto 1.0 3.0 100] + set result [::math::statistics::estimate-pareto $rnumbers] +} -result {1.000519 3.668162 0.3668162} test "kruskal-wallis-1.0" "Test analysis Kruskal-Wallis" -match tolerant -body { ::math::statistics::analyse-Kruskal-Wallis {6.4 6.8 7.2 8.3 8.4 9.1 9.4 9.7} {2.5 3.7 4.9 5.4 5.9 8.1 8.2} {1.3 4.1 4.9 5.2 5.5 8.2} } -result {9.83627087199 0.00731275323967} test "kruskal-wallis-1.1" "Test test Kruskal-Wallis" -match tolerant -body { Index: modules/math/wilcoxon.tcl ================================================================== --- modules/math/wilcoxon.tcl +++ modules/math/wilcoxon.tcl @@ -1,11 +1,10 @@ # statistics_new.tcl -- # Implementation of the Wilcoxon test: test if the medians # of two samples are the same # -package require math::statistics # test-Wilcoxon # Compute the statistic that indicates if the medians of two # samples are the same #