Tcl Library Source Code

2015-04-29  Arjen Markus <[email protected]>
	* 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 <[email protected]>
	* statistics.tcl: Bump version to 1.0 - Aku found the cause of earlier problems

2015-04-26  Arjen Markus <[email protected]>
	* 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 <[email protected]>
	* 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)
	* special.test: Adding test case for this new function
	* special.man: Describing invnorm plus a correction in the overview (ierfc_n is not implemented)

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


dd. 26 october 2005, Arjen Markus

# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
# 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.1

# math::calculus --
#    Namespace for the commands

namespace eval ::math::calculus {

    namespace import ::math::interpolate::neville

           0.1047900103222502e+00    0.1406532597155259e+00
           0.1690047266392679e+00    0.1903505780647854e+00
           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

    #

    }

    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 [expr {pow((0.2e+3*$abserr/$resasc),1.5e+00)}]]}]
    }
    if { $resabs > $uflow/(0.5e+02*$epmach) } {
        set abserr [Max [expr {($epmach*0.5e+02)*$resabs}] $abserr]
    }

    return [list $result $abserr $resabs $resasc]
}

# qk15 --
#     Apply the QK15 rule to an interval and return the estimated integral

        set abserr 0.0
        set resabs 0.0
        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)}]

            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]
}

    if { $m < $n } {
        set U {}
        incr m -1
        foreach row $A {
            lappend U [lrange $row 0 $m]
        }
        #puts $U
    }
    return [list $U $S $V]
}

# eigenvectorsSVD --
#     Determine the eigenvectors and eigenvalues of a real
#     symmetric matrix via the SVD

# mvreglin.tcl --
#     Addition to the statistics package
#     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


# ::math::statistics --
#     This file adds:
#     mvlinreg = Multivariate Linear Regression
#
namespace eval ::math::statistics {
    variable epsilon 1.0e-7

#

# ::math::statistics --
#   Namespace holding the procedures and variables
#
namespace eval ::math::statistics {

    namespace export pdf-normal pdf-uniform pdf-lognormal \
	    pdf-exponential \
	    cdf-normal cdf-uniform cdf-lognormal \
	    cdf-exponential \
	    cdf-students-t \
	    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 \
	    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
    variable twopi               [expr {2.0*acos(-1.0)}]
    variable pi                  [expr {acos(-1.0)}]

    if { $stdev <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $NEGSTDEV $NEGSTDEV
    }

    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 --
#    Return the probabilities belonging to a uniform distribution
#    (parameters as minimum/maximum)

	return -code error -errorcode ARG -errorinfo $NEGSTDEV $NEGSTDEV
    }

    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
}

        }
    }

    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:
#    m         First number of degrees of freedom
#    n         Second number of degrees of freedom

    if {$b <= 0.0} {
        return -code error "Value out of range in Beta density: b = $b, must be > 0"
    }
    #
    # Corner cases ... need to check these!
    #
    if {$x == 0.0} {
        return [expr {$a > 1.0? 0.0 : Inf}]
    }
    if {$x == 1.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)}]

    set term [expr {$term1 + $term2}]
    if { $term > -200.0 } {

    } elseif {$retval > 1.0} {
        set retval 1.0
    }

    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:
#    alpha     Shape parameter
#    beta      Rate parameter

# Note:
#    Implemented by Eric Kemp-Benedict, 2008
#
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
#    number    Number of values to return
#
# Result:
#    List of random values
#
# Note:
#    Implemented by Eric Kemp-Benedict, 2007
#    Generate a list of gamma-distributed random deviates

        }

        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]) } {

    #

package ifneeded math::fourier           1.0.2 [list source [file join $dir fourier.tcl]]

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.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.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]]

if {![package vsatisfies [package provide Tcl] 8.5]} {return}

[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::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]
[arg_def list data1] - First list of data
[arg_def list data2] - Second list of data
[list_end]
[para]

[call [cmd ::math::statistics::interval-mean-stdev] [arg data] [arg confidence]]
Return the interval containing the mean value and one
containing the standard deviation with a certain
level of confidence (assuming a normal distribution)

[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 alpha]]
Test whether the mean value of a sample is in accordance with the
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 alpha] - Probability level (0.95 or 0.99 for instance)
[list_end]
[para]


[call [cmd ::math::statistics::test-normal] [arg data] [arg significance]]
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 data] - List of raw data values
[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
points.

[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
[para]
[arg_def float confidence] - Confidence level (0.95 or 0.99 for instance) or a list of confidence levels.

}]

[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 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:
[list_begin itemized]
[item]

[list_begin definitions]

[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]
[para]

distribution with given degrees of freedom

[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] - 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]
[para]

[arg_def int degrees] - Number of degrees of freedom
[arg_def float value] - Value for which the probability is required
[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.

[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
[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).

[list_begin arguments]
[arg_def float mu] - Mean number of occurrences
[arg_def int k] - Number of occurences
[list_end]
[para]

[call [cmd ::math::statistics::cdf-beta] [arg a] [arg b] [arg value]]
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] - 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]
[para]

[arg_def float xmax] - Maximum value of the distribution
[arg_def int number] - Number of values to be returned
[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.

[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
[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.

[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.

[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.

[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.

[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]
[arg_def float xmin] - Minimum value of the distribution
[arg_def float xmax] - Maximum value of the distribution

[list_begin arguments]
[arg_def float a] - First shape parameter
[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"]
The data manipulation procedures act on lists or lists of lists:

# version 0.5:   added the population standard deviation and variance,
#                as suggested by Dimitrios Zachariadis
# 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 1.0
package require math

if {![llength [info commands ::lrepeat]]} {
    # Forward portability, emulate lrepeat
    proc ::lrepeat {n args} {
	if {$n < 1} {
	    return -code error "must have a count of at least 1"

#

namespace eval ::math::statistics {
    #
    # Safer: change to short procedures
    #
    namespace export mean min max number var stdev pvar pstdev basic-stats corr \
	    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-Duckworth
    #
    # Error messages
    #
    variable NEGSTDEV   {Zero or negative standard deviation}
    variable TOOFEWDATA {Too few or invalid data}
    variable OUTOFRANGE {Argument out of range}

	    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
	}

	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
	}

}

# test-normal --
#     Test for normality (using method Lilliefors)
#
# Arguments:
#     data            Values that need to be tested
#     significance    Level at which the discrepancy from normality is tested
#
# Result:
#     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 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($significance-0.20) < 0.0001 } {
        set Dcrit 0.741
    }
    if { abs($significance-0.15) < 0.0001 } {
        set Dcrit 0.775
    }
    if { abs($significance-0.10) < 0.0001 } {
        set Dcrit 0.819
    }
    if { abs($significance-0.05) < 0.0001 } {
        set Dcrit 0.895
    }
    if { abs($significance-0.01) < 0.0001 } {
        set Dcrit 1.035
    }
    if { $Dcrit != "--" } {
        return [expr {$D > $Dcrit ? 1 : 0 }]
    } else {
        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 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
#    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 alpha } {
    variable NEGSTDEV
    variable TOOFEWDATA

    if { $est_stdev <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $NEGSTDEV $NEGSTDEV
    }

    set allstats        [BasicStats all $data]

    set alpha2          [expr {(1.0+$alpha)/2.0}]

    set sample_mean     [lindex $allstats 0]
    set sample_number   [lindex $allstats 3]

    if { $sample_number > 1 } {
	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<$alpha2}]

    } else {
	return -code error -errorcode DATA -errorinfo $TOOFEWDATA $TOOFEWDATA
    }
}

# interval-mean-stdev --

        if { $range < $rlower } { lappend result $i }
        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]
source [file join [file dirname [info script]] plotstat.tcl]
source [file join [file dirname [info script]] liststat.tcl]

} -result {1 1 8 0}

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
#
test "Quantiles-1.0" "Quantiles - raw data" -match tolerant -body {

#
test "Testnormal-1.0" "Determine normality statistic for birth weight data" -match tolerant -body {
    ::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 - 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.20
} -result 1

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.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
test "Testmultivar-1.0" "Ordinary multivariate regression - three independent variables" \
        -match tolerant -body {

        [::math::statistics::cdf-normal  0   1   1] \
        [::math::statistics::cdf-normal  0.0 1.0 1.0] \
        [::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]]
} -result {0.00263194027271168 0.0302770403110644 2.62677891379834e-07}

        [::math::statistics::cdf-beta 1000 1000 0.7] \
        [::math::statistics::cdf-beta 2 3 0.6]]
} -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
#
# 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

test "random-numbers-1.3" "Test random-gamma" -body {

        if { $r < 0.0 || $r > 1.0 } {
            result 0
            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 {
    ::math::statistics::test-Kruskal-Wallis 0.95 {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 1

# statistics_new.tcl --
#     Implementation of the Wilcoxon test: test if the medians
#     of two samples are the same
#



# test-Wilcoxon
#     Compute the statistic that indicates if the medians of two
#     samples are the same
#
# Arguments:
#     sample_a       List of values in the first sample