cdf-exponential cdf-triangular cdf-symmetric-triangular \
	    cdf-students-t \
	    random-normal random-uniform random-lognormal \
	    random-exponential random-triangular \
	    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-F \
	    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 bootstrap estimate-exponential


    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 { $x < 0.0 } { return 0.0 }
    if { $x > 700.0*$mean } { return 0.0 }

    set prob [expr {exp(-$x/double($mean))/$mean}]























































































    return $prob
}


# cdf-normal --
#    Return the cumulative probability belonging to a normal distribution
#
................................................................................
    if { $x > 30.0*$mean } { return 1.0 }

    set prob [expr {1.0-exp(-$x/double($mean))}]

    return $prob
}































































































# Inverse-cdf-uniform --
#    Return the argument belonging to the cumulative probability
#    for a uniform distribution (parameters as minimum/maximum)
#
# Arguments:
#    pmin      Minimum of the distribution
................................................................................
    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
#
................................................................................
    }

    set parameter [expr {$sum/double($count)}]
    set stdev     [expr {$parameter / sqrt($count)}]

    return [list $parameter $stdev]
}

































































# empirical-distribution --
#    Determine the empirical distribution
#
# Arguments:
#    values    Values that are to be examined
#
................................................................................
    foreach z    {4.417 3.891 3.291 2.576 2.241 1.960 1.645 1.150 0.674
    0.319 0.126 0.063 0.0125} \
	    pexp {1.e-5 1.e-4 1.e-3 1.e-2 0.025 0.050 0.100 0.250 0.500
    0.750 0.900 0.950 0.990 } {
	set prob [::math::statistics::Cdf-toms322 1 5000 [expr {$z*$z}]]
	puts "$z - $pexp - [expr {1.0-$prob}]"
    }

    puts "Normal distribution (inverted; one-tailed)"
    foreach p {0.001 0.01 0.1 0.25 0.5 0.75 0.9 0.99 0.999} {
	puts "$p - [::math::statistics::Inverse-cdf-normal 0.0 1.0 $p]"
    }
    puts "Normal random variables"
    set rndvars [::math::statistics::random-normal 1.0 2.0 20]
    puts $rndvars

	    cdf-exponential cdf-triangular cdf-symmetric-triangular \
	    cdf-students-t \
	    random-normal random-uniform random-lognormal \
	    random-exponential random-triangular \
	    histogram-uniform \
	    pdf-gamma pdf-poisson pdf-chisquare pdf-students-t pdf-beta \
	    pdf-weibull pdf-gumbel pdf-pareto pdf-cauchy \
	    pdf-laplace pdf-kumaraswamy pdf-negative-binomial \
	    cdf-gamma cdf-poisson cdf-chisquare cdf-beta cdf-F \
	    cdf-weibull cdf-gumbel cdf-pareto cdf-cauchy \
	    cdf-laplace cdf-kumaraswamy cdf-negative-binomial \
	    random-gamma random-poisson random-chisquare random-students-t random-beta \
	    random-weibull random-gumbel random-pareto random-cauchy \
	    random-laplace random-kumaraswamy random-negative-binomial \
	    incompleteGamma incompleteBeta \
	    estimate-pareto empirical-distribution bootstrap estimate-exponential \
	    estimate-laplace estimate-negative-binomial

    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 { $x < 0.0 } { return 0.0 }
    if { $x > 700.0*$mean } { return 0.0 }

    set prob [expr {exp(-$x/double($mean))/$mean}]

    return $prob
}


# pdf-laplace --
#    Return the probabilities belonging to a Laplace
#    distribution
#
# Arguments:
#    mean     Mean of the distribution
#    scale    Scale (the spreading) of the distribution
#    x        Value for which the probability must be determined
#
# Result:
#    Probability of value x under the given distribution
#
proc ::math::statistics::pdf-laplace { mean scale x } {
    variable NEGSTDEV
    variable OUTOFRANGE

    if { $scale <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: scale must be positive"
    }

    set prob [expr {exp(-($x-$mean)/double($scale))/(2.0*$scale)}]

    return $prob
}


# pdf-kumaraswamy --
#    Return the probabilities belonging to a Kumaraswamy
#    distribution (akin to the Beta distribution, but tractable)
#
#    Arguments:
#    a         First parameter of the Kumaraswamy distribution
#    b         Second parameter of the Kumaraswamy distribution
#    x         Value of variate
#
# Result:
#    Probability of value x under the given distribution
#
proc ::math::statistics::pdf-kumaraswamy { a b x } {
    variable OUTOFRANGE

    if { $a <= 0.0 || $b <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameters a and b must be positive"
    }

    set prob [expr {$a * $b * $x**($a-1) * (1.0 -$x**$a)**($b-1)}]

    return $prob
}


# pdf-negative-binomial --
#    Return the probability belonging to a negative binomial
#    distribution
#
#    Arguments:
#    r         Allowed number of failures for the distribution
#    p         Probability of success for the negative bionomial distribution
#    k         Value of variate (integer)
#
# Result:
#    Probability of k successes under the given distribution
#
proc ::math::statistics::pdf-negative-binomial { r p k } {
    variable OUTOFRANGE
    variable INVALID

    if { $p < 0.0 || $p >= 1.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameter p must be non-negative and lower than 1"
    }

    if { int($r) != $r || $r < 1 } {
	return -code error -errorcode ARG -errorinfo $INVALIDE \
		"$INVALID: parameter r must be a positive integer"
    }

    set coeff [::math::choose [expr {$k+$r-1}] $k]
    set prob  [expr {$coeff * (1.0 - $p)**$r * $p ** $k}]

    return $prob
}


# cdf-normal --
#    Return the cumulative probability belonging to a normal distribution
#
................................................................................
    if { $x > 30.0*$mean } { return 1.0 }

    set prob [expr {1.0-exp(-$x/double($mean))}]

    return $prob
}


# cdf-laplace --
#    Return the cumulative probabilities belonging to a Laplace
#    distribution
#
# Arguments:
#    mean     Mean of the distribution
#    scale    Scale (the spreading) of the distribution
#    x        Value for which the probability must be determined
#
# Result:
#    Cumulative probability of value x under the given distribution
#
proc ::math::statistics::cdf-laplace { mean scale x } {
    variable NEGSTDEV
    variable OUTOFRANGE

    if { $scale <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: scale must be positive"
    }

    if { $x < $mean } {
        set prob [expr {0.5 * exp(($x-$mean)/double($scale))}]
    } else {
        set prob [expr {1.0 - 0.5 * exp(($mean-$x)/double($scale))}]
    }

    return $prob
}


# cdf-kumaraswamy --
#    Return the cumulative probabilities belonging to a Kumaraswamy
#    distribution (akin to the Beta distribution, but tractable)
#
#    Arguments:
#    a         First parameter of the Kumaraswamy distribution
#    b         Second parameter of the Kumaraswamy distribution
#    x         Value of variate
#
# Result:
#    Cumulative probability of value x under the given distribution
#
proc ::math::statistics::cdf-kumaraswamy { a b x } {
    variable OUTOFRANGE

    if { $a <= 0.0 || $b <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameters a and b must be positive"
    }

    set prob [expr {1.0 - (1.0-$x**$a) ** $b}]

    return $prob
}


# cdf-negative-binomial --
#    Return the cumulative probability for a negative binomial distribution
#
#    Arguments:
#    r         Allowed number of failures for the distribution
#    p         Probability of success for the negative bionomial distribution
#    k         Value of variate (integer)
#
# Result:
#    Cumulative probability for up to k successes under the given distribution
#
proc ::math::statistics::cdf-negative-binomial { r p k } {
    variable OUTOFRANGE
    variable INVALID

    if { $p < 0.0 || $p >= 1.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameter p must be non-negative and lower than 1"
    }

    if { int($r) != $r || $r < 1 } {
	return -code error -errorcode ARG -errorinfo $INVALIDE \
		"$INVALID: parameter r must be a positive integer"
    }

    set sum 0.0

    for { set i 0 } { $i <= $k } { incr i } {
        set prob [pdf-negative-binomial $r $p $i]
        set sum  [expr {$sum + $prob}]
    }

    return $sum
}


# Inverse-cdf-uniform --
#    Return the argument belonging to the cumulative probability
#    for a uniform distribution (parameters as minimum/maximum)
#
# Arguments:
#    pmin      Minimum of the distribution
................................................................................
    set retval {}
    for {set i 0} {$i < $number} {incr i} {
        lappend retval [expr {$location + $scale * tan( $pi * (rand() - 0.5))}]
    }
    return $retval
}


# random-laplace --
#    Generate a list of Laplace distributed deviates
#
# Arguments:
#    mean     Mean of the distribution
#    scale    Scale (the spreading) of the distribution
#    number   Number of values to return
#
# Result:
#    List of random numbers
#
proc ::math::statistics::random-laplace { mean scale number } {
    variable NEGSTDEV
    variable OUTOFRANGE

    if { $scale <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: scale must be positive"
    }

    set retval {}
    for {set i 0} {$i < $number} {incr i} {
        set p [expr {rand()}]
        if { $p < 0.5 } {
            set x [expr {$mean + $scale * log(1.0 - 2.0*abs($p-0.5))}]
        } else {
            set x [expr {$mean - $scale * log(1.0 - 2.0*abs($p-0.5))}]
        }
        lappend retval $x
    }

    return $retval
}


# random-kumaraswamy --
#    Generate a list of Kumaraswamy distributed deviates
#
#    Arguments:
#    a         First parameter of the Kumaraswamy distribution
#    b         Second parameter of the Kumaraswamy distribution
#    number    Number of values to return
#
# Result:
#    List of random numbers
#
proc ::math::statistics::random-kumaraswamy { a b number } {
    variable OUTOFRANGE

    if { $a <= 0.0 || $b <= 0.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameters a and b must be positive"
    }

    set ra [expr {1.0 / $a}]
    set rb [expr {1.0 / $b}]

    set retval {}
    for {set i 0} {$i < $number} {incr i} {
        set p [expr {rand()}]
        set x [expr {( 1.0 - (1.0-$p) ** $rb ) ** $ra}]

        lappend retval $x
    }

    return $retval
}


# random-negative-binomial --
#    Generate a list of deviates according to the negative binomial distribution
#
#    Arguments:
#    r         Allowed number of failures for the distribution
#    p         Probability of success for the negative bionomial distribution
#    number    Number of values to return
#
# Result:
#    List of random numbers
#
proc ::math::statistics::random-negative-binomial { r p number } {
    variable OUTOFRANGE
    variable INVALID

    if { $p < 0.0 || $p >= 1.0 } {
	return -code error -errorcode ARG -errorinfo $OUTOFRANGE \
		"$OUTOFRANGE: parameter p must be non-negative and lower than 1"
    }

    if { int($r) != $r }} $r < 1} {
	return -code error -errorcode ARG -errorinfo $INVALIDE \
		"$INVALID: parameter r must be a positive integer"
    }

    set retval {}
    for {set i 0} {$i < $number} {incr i} {
        set success 0
        set failure 0

        while { $failure < $r } {
            if { rand() <= $p } {
                incr success
            } else {
                incr failure
            }
        }

        lappend retval $success
    }

    return $retval
}


# estimate-pareto --
#    Estimate the parameters of a Pareto distribution
#
# Arguments:
#    values    Values that are supposed to be distributed according to Pareto
#
................................................................................
    }

    set parameter [expr {$sum/double($count)}]
    set stdev     [expr {$parameter / sqrt($count)}]

    return [list $parameter $stdev]
}


# estimate-laplace --
#    Estimate the parameter of a Laplace distribution
#
# Arguments:
#    values    Values that are supposed to be Laplace distributed
#
# Result:
#    Estimates of respectively the mean and the scale of the Laplace distribution
#    (See https://en.wikipedia.org/wiki/Laplace_distribution)
#
# Note:
#    According to Wikipedia the estimators are maximum-likelihood estimators
#
proc ::math::statistics::estimate-laplace { values } {

    set mean [median $values]

    set sum   0.0
    set count 0

    foreach v $values {
        if { $v != "" } {
            set  sum [expr {$sum + abs($v-$mean)}]
            incr count
        }
    }

    set scale [expr {$sum/double($count)}]

    return [list $mean $scale]
}


# estimate-negative-binomial --
#    Estimate the parameter p of a negative binomial distribution,
#    given the allowed number of failures
#
# Arguments:
#    r         Allowed number of failures
#    values    Values that are supposed to be distributed according to a negative binomial distribution
#
# Result:
#    Estimate of the probability of success for the distribution
#
# Note:
#    According to Wikipedia the estimators are maximum-likelihood estimators
#
proc ::math::statistics::estimate-negative-binomial { r values } {

    set sum   0.0
    set count 0

    foreach v $values {
        if { $v != "" } {
            set  sum [expr {$sum + $v}]
            incr count
        }
    }

    return [expr {$sum/double($count * $r + $sum)}]
}


# empirical-distribution --
#    Determine the empirical distribution
#
# Arguments:
#    values    Values that are to be examined
#
................................................................................
    foreach z    {4.417 3.891 3.291 2.576 2.241 1.960 1.645 1.150 0.674
    0.319 0.126 0.063 0.0125} \
	    pexp {1.e-5 1.e-4 1.e-3 1.e-2 0.025 0.050 0.100 0.250 0.500
    0.750 0.900 0.950 0.990 } {
	set prob [::math::statistics::Cdf-toms322 1 5000 [expr {$z*$z}]]
	puts "$z - $pexp - [expr {1.0-$prob}]"
    }

    puts "Normal distribution (inverted; one-tailed)"
    foreach p {0.001 0.01 0.1 0.25 0.5 0.75 0.9 0.99 0.999} {
	puts "$p - [::math::statistics::Inverse-cdf-normal 0.0 1.0 $p]"
    }
    puts "Normal random variables"
    set rndvars [::math::statistics::random-normal 1.0 2.0 20]
    puts $rndvars

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}
package ifneeded math::calculus          0.8.1 [list source [file join $dir calculus.tcl]]
# statistics depends on linearalgebra (for multi-variate linear regression).
# statistics depends on optimize (for logistic regression).
package ifneeded math::statistics        1.3.1 [list source [file join $dir statistics.tcl]]
package ifneeded math::linearalgebra     1.1.6 [list source [file join $dir linalg.tcl]]
package ifneeded math::calculus::symdiff 1.0.1 [list source [file join $dir symdiff.tcl]]
package ifneeded math::bigfloat          2.0.2 [list source [file join $dir bigfloat2.tcl]]
package ifneeded math::numtheory         1.1.1 [list source [file join $dir numtheory.tcl]]
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.3.1 [list source [file join $dir geometry.tcl]]
package ifneeded math::trig              1.0   [list source [file join $dir trig.tcl]]
package ifneeded math::quasirandom       1.0   [list source [file join $dir quasirandom.tcl]]
package ifneeded math::special           0.4.0 [list source [file join $dir special.tcl]]

if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0.1 [list source [file join $dir exact.tcl]]
package ifneeded math::PCA               1.0   [list source [file join $dir pca.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}
package ifneeded math::calculus          0.8.1 [list source [file join $dir calculus.tcl]]
# statistics depends on linearalgebra (for multi-variate linear regression).
# statistics depends on optimize (for logistic regression).
package ifneeded math::statistics        1.5.0 [list source [file join $dir statistics.tcl]]
package ifneeded math::linearalgebra     1.1.6 [list source [file join $dir linalg.tcl]]
package ifneeded math::calculus::symdiff 1.0.1 [list source [file join $dir symdiff.tcl]]
package ifneeded math::bigfloat          2.0.2 [list source [file join $dir bigfloat2.tcl]]
package ifneeded math::numtheory         1.1.1 [list source [file join $dir numtheory.tcl]]
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.3.1 [list source [file join $dir geometry.tcl]]
package ifneeded math::trig              1.0   [list source [file join $dir trig.tcl]]
package ifneeded math::quasirandom       1.0   [list source [file join $dir quasirandom.tcl]]
package ifneeded math::special           0.4.0 [list source [file join $dir special.tcl]]

if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0.1 [list source [file join $dir exact.tcl]]
package ifneeded math::PCA               1.0   [list source [file join $dir pca.tcl]]

# version 0.9.3: added histogram-alt, corrected test-normal
# version 1.0:   added test-anova-F
# version 1.0.1: correction in pdf-lognormal and cdf-lognormal
# version 1.1:   added test-Tukey-range and test-Dunnett
# version 1.3:   added wasserstein-distance, kl-divergence and logit regression

package require Tcl 8.5 ; # 8.5+ feature in test-anova-F and others: **-operator
package provide math::statistics 1.3.1
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"
................................................................................
	    test-Duckworth test-anova-F test-Tukey-range test-Dunnett
    #
    # Error messages
    #
    variable NEGSTDEV   {Zero or negative standard deviation}
    variable TOOFEWDATA {Too few or invalid data}
    variable OUTOFRANGE {Argument out of range}


    #
    # Coefficients involved
    #
    variable factorNormalPdf
    set factorNormalPdf [expr {sqrt(8.0*atan(1.0))}]

# version 0.9.3: added histogram-alt, corrected test-normal
# version 1.0:   added test-anova-F
# version 1.0.1: correction in pdf-lognormal and cdf-lognormal
# version 1.1:   added test-Tukey-range and test-Dunnett
# version 1.3:   added wasserstein-distance, kl-divergence and logit regression

package require Tcl 8.5 ; # 8.5+ feature in test-anova-F and others: **-operator
package provide math::statistics 1.5.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"
................................................................................
	    test-Duckworth test-anova-F test-Tukey-range test-Dunnett
    #
    # Error messages
    #
    variable NEGSTDEV   {Zero or negative standard deviation}
    variable TOOFEWDATA {Too few or invalid data}
    variable OUTOFRANGE {Argument out of range}
    variable INVALID    {Argument invalid}

    #
    # Coefficients involved
    #
    variable factorNormalPdf
    set factorNormalPdf [expr {sqrt(8.0*atan(1.0))}]

# 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:
................................................................................
#
# Result:
#     Rank correlation coefficient
#
proc ::math::statistics::spearman-rank {sample_a sample_b} {
    return [lindex [spearman-rank-extended $sample_a $sample_b] 0]
}

# wilcoxon.tcl --
#     Implementation of the Wilcoxon test: test if the medians
#     of two samples are the same
#
#     Also: Levene's and Brown-Forsythe's test
#

# test-Wilcoxon
#     Compute the statistic that indicates if the medians of two
#     samples are the same
#
# Arguments:
................................................................................
#
# Result:
#     Rank correlation coefficient
#
proc ::math::statistics::spearman-rank {sample_a sample_b} {
    return [lindex [spearman-rank-extended $sample_a $sample_b] 0]
}

# test-Levene --
#     Compute the Levene statistic that indicates if the variances of
#     groups of data are the same
#
# Arguments:
#     groups         List of groups of values to be examined
#
# Result:
#     Statistic for the test (an F statistic with k-1, N-k degrees
#     of freedom - k the number of groups and N the total number
#     of values)
#     The test uses the mean of the values in the groups.
#
proc ::math::statistics::test-Levene {groups} {

    return [Test-Levene-Brown-Forsythe 0 $groups]
}

# test-Brown-Forsythe --
#     Compute the Brown-Forsythe statistic that indicates if the variances of
#     groups of data are the same
#
# Arguments:
#     groups         List of groups of values to be examined
#
# Result:
#     Statistic for the test (an F statistic with k-1, N-k degrees
#     of freedom - k the number of groups and N the total number
#     of values)
#     The test uses the median of the values in the groups.
#
proc ::math::statistics::test-Brown-Forsythe {groups} {

    return [Test-Levene-Brown-Forsythe 1 $groups]
}

# Test-Levene-Brown-Forsythe --
#     Compute either the Levene or the Brown-Forsythe statistic that indicates
#     if the variances of groups of data are the same
#
# Arguments:
#     choice         Which of the two versions
#     groups         List of groups of values to be examined
#
# Result:
#     Statistic for the test
#     The test uses either the mean or the median of the values in the groups.
#
proc ::math::statistics::Test-Levene-Brown-Forsythe {choice groups} {

    #
    # Compute the deviations from the mean/median within each group
    #
    set alldevs {}
    set zscores {}
    set zmeans  {}
    foreach group $groups {
        if { $choice } {
            set zm [median $group]
        } else {
            set zm [mean $group]
        }
        set zgroup {}
        foreach element $group {
            lappend zgroup [expr {abs($element-$zm)}]
        }

        set alldevs [concat $alldevs $zgroup]
        lappend zscores $zgroup
        lappend zmeans  [mean $zgroup]
    }

    set zoverall [mean $alldevs]

    set ndata   [llength $alldevs]
    set ngroups [llength $groups]

    #
    # Compute the numerator of the statistic
    #
    set sumsqmeans 0.0

    foreach zm $zmeans group $groups {
        set n          [llength $group]
        set sumsqmeans [expr { $sumsqmeans + $n * ($zm - $zoverall)**2 }]
    }

    #
    # Compute the denominator
    #
    set sumsqpergroup 0.0

    foreach zm $zmeans zs $zscores {
        set sumsq 0.0
        foreach z $zs {
            set sumsq [expr {$sumsq + ($z-$zm)**2}]
        }

        set sumsqpergroup [expr { $sumsqpergroup + $sumsq }]
    }

    #
    # Finally, the statistic
    #

    return [expr { ($ndata-$ngroups) * $sumsqmeans / double( ($ngroups-1) * $sumsqpergroup ) }]
}