Tcl Library Source Code

2014-01-18  Arjen Markus <[email protected]>
	* statistics.tcl: Added stat_kernel.tcl
	* stat_kernel.tcl: Implements a straightforward kernel density estimation procedure
	* statistics.man: Describe the kernel denstity estimation procedure, moved the description of several
	                  tests to the general section
	* pkgIndex.tcl: Bumped version of statistics package to 0.9

2013-12-20  Arjen Markus <[email protected]>
	* interpolate.tcl: [Ticket 843c2257d2] Added special case for points coincident with the data points
	* interpolate.test: [Ticket 843c2257d2] Added test case for coincident points

2013-12-17  Andreas Kupries  <[email protected]>

	* decimal.man: Fixed missing requirement of the package

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   [list source [file join $dir statistics.tcl]]
package ifneeded math::optimize          1.0   [list source [file join $dir optimize.tcl]]
package ifneeded math::calculus          0.7.2 [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.4 [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]]

# stat_kernel.tcl --
#
#    Part of the statistics package for basic statistical analysis
#    Based on http://en.wikipedia.org/wiki/Kernel_(statistics) and
#             http://en.wikipedia.org/wiki/Kernel_density_estimation
#
# version 0.1:   initial implementation, january 2014

# kernel-density --
#     Estimate the probability density using the kernel density
#     estimation method
#
# Arguments:
#     data            List of univariate data
#     args            List of options in the form of keyword-value pairs:
#                     -weights weights: per data point the weight
#                     -bandwidth value: bandwidth to be used for the estimation
#                     -number value: number of bins to be returned
#                     -interval {begin end}: begin and end of the interval for
#                         which the density is returned
#                     -kernel function: kernel to be used (gaussian, cosine,
#                         epanechnikov, uniform, triangular, biweight,
#                         logistic)
#                     For all options more or less sensible defaults are
#                     provided.
#
# Result:
#     A list of the bin centres, a list of the corresponding density
#     estimates and a list containing several computational parameters:
#     begin and end of the interval, mean, standard deviation and bandwidth
#
# Note:
#     The conditions for the kernel function are fairly weak:
#     - It should integrate to 1
#     - It should be symmetric around 0
#
#     As for the implementation in Tcl: it should be reachable in the
#     ::math::statistics namespace. As a consequence, you can define
#     your own kernel function too. Hence there is no check.
#
proc ::math::statistics::kernel-density {data args} {

    #
    # Determine the basic statistics
    #
    set basicStats [BasicStats all $data]

    set mean       [lindex $basicStats 0]
    set ndata      [lindex $basicStats 3]
    set stdev      [lindex $basicStats 4]

    if { $ndata < 1 } {
        return -code error -errorcode ARG -errorinfo "Too few actual data"
    }

    #
    # Get the options (providing defaults as needed)
    #
    set opt(-weights)   {}
    set opt(-number)    100
    set opt(-kernel)    gaussian

    #
    # The default bandwidth is set via a simple expression, which
    # is supposed to be optimal for the Gaussian kernel.
    # Perhaps a more sophisticated method should be provided as well
    #
    set opt(-bandwidth) [expr {1.06 * $stdev / pow($ndata,0.2)}]

    #
    # The default interval is derived from the mean and the
    # standard deviation
    #
    set opt(-interval) [list [expr {$mean - 3.0 * $stdev}] [expr {$mean + 3.0 * $stdev}]]

    #
    # Retrieve the given options from $args
    #
    if { [llength $args] % 2 != 0 } {
        return -code error -errorcode ARG -errorinfo "The options must all have a value"
    }
    array set opt $args

    #
    # Elementary checks
    #
    if { $opt(-bandwidth) <= 0.0 } {
        return -code error -errorcode ARG -errorinfo "The bandwidth must be positive: $opt(-bandwidth)"
    }

    if { $opt(-number) <= 0.0 } {
        return -code error -errorcode ARG -errorinfo "The number of bins must be positive: $opt(-number)"
    }

    if { [llength [info proc $opt(-kernel)]] == 0 } {
        return -code error -errorcode ARG -errorinfo "Unknown kernel function: $opt(-kernel)"
    }

    #
    # Construct the weights
    #
    if { [llength $opt(-weights)] > 0 } {
        if { [llength $data] != [llength $opt(-weights)] } {
            return -code error -errorcode ARG -errorinfo "The list of weights must match the data"
        }

        set sum 0.0
        foreach d $data w $opt(-weights) {
            if { $d != {} } {
                set sum [expr {$sum + $w}]
            }
        }
        set scale [expr {1.0/$sum/$ndata}]

        set weight {}
        foreach w $opt(-weights) {
            if { $d != {} } {
                lappend weight [expr {$w / $scale}]
            } else {
                lappend weight {}
            }
        }
    } else {
        set weight [lrepeat $ndata [expr {1.0/$ndata}]]
    }

    #
    # Construct the centres of the bins
    #
    set xbegin [lindex $opt(-interval) 0]
    set xend   [lindex $opt(-interval) 1]
    set dx     [expr {($xend - $xbegin) / $opt(-number)}]
    set xb     [expr {$xbegin + 0.5 * $dx}]
    set xvalue {}
    for {set i 0} {$i < $opt(-number)} {incr i} {
        lappend xvalue [expr {$xb + $i * $dx}]
    }

    #
    # Construct the density function
    #
    set density {}
    set scale   [expr {$opt(-bandwidth)}]
    foreach x $xvalue {
        set sum 0.0
        foreach d $data w $weight {
            if { $d != {} } {
                set kvalue [$opt(-kernel) [expr {$scale * ($x-$d)}]]
                set sum [expr {$sum + $w * $kvalue}]
            }
        }
        lappend density [expr {$sum * $scale}]
    }

    #
    # Return the result
    #
    return [list $xvalue $density [list $xbegin $xend $mean $stdev $scale]]
}

# gaussian, uniform, triangular, epanechnikov, biweight, cosine, logistic --
#    The Gaussian kernel
#
# Arguments:
#    x            (Scaled) argument
#
# Result:
#    Value of the kernel
#
# Note:
#    The standard deviation is 1.
#
proc ::math::statistics::gaussian {x} {
    return [expr {exp(-0.5*$x*$x) / sqrt(2.0*acos(-1.0))}]
}
proc ::math::statistics::uniform {x} {
    if { abs($x) < 1.0 } {
        return 0.5
    } else {
        return 0.0
    }
}
proc ::math::statistics::triangular {x} {
    if { abs($x) < 1.0 } {
        return [expr {1.0 - abs($x)}]
    } else {
        return 0.0
    }
}
proc ::math::statistics::epanechnikov {x} {
    if { abs($x) < 1.0 } {
        return [expr {0.75 * (1.0 - abs($x)*abs($x))}]
    } else {
        return 0.0
    }
}
proc ::math::statistics::biweight {x} {
    if { abs($x) < 1.0 } {
        return [expr {0.9375 * pow((1.0 - abs($x)*abs($x)),2)}]
    } else {
        return 0.0
    }
}
proc ::math::statistics::cosine {x} {
    if { abs($x) < 1.0 } {
        return [expr {0.25 * acos(-1.0) * cos(0.5 * acos(-1.0) * $x)}]
    } else {
        return 0.0
    }
}
proc ::math::statistics::logistic {x} {
    return [expr {1.0 / (exp($x) + 2.0 + exp(-$x))}]
}

[manpage_begin math::statistics n 0.9]
[keywords {data analysis}]
[keywords mathematics]
[keywords statistics]
[moddesc {Tcl Math Library}]
[titledesc {Basic statistical functions and procedures}]
[category  Mathematics]
[require Tcl 8.4]
[require math::statistics 0.9]
[description]
[para]

The [package math::statistics] package contains functions and procedures for
basic statistical data analysis, such as:

[list_begin itemized]

Returns a list of subsamples (their indices) that indeed violate the
limits.
[list_begin arguments]
[arg_def list control] - Control limits as returned by the "control-Rchart" procedure
[arg_def list data] - List of observed data
[list_end]
[para]

[call [cmd ::math::statistics::test-Kruskal-Wallis] [arg confidence] [arg args]]
Check if the population medians of two or more groups are equal with a
given confidence level, using the Kruskal-Wallis test.

[list_begin arguments]
[arg_def float confidence] - Confidence level to be used (0-1)
[arg_def list args] - Two or more lists of data
[list_end]
[para]

[call [cmd ::math::statistics::analyse-Kruskal-Wallis] [arg args]]
Compute the statistical parameters for the Kruskal-Wallis test.
Returns the Kruskal-Wallis statistic and the probability that that
value would occur assuming the medians of the populations are
equal.

[list_begin arguments]
[arg_def list args] - Two or more lists of data
[list_end]
[para]

[call [cmd ::math::statistics::group-rank] [arg args]]
Rank the groups of data with respect to the complete set.
Returns a list consisting of the group ID, the value and the rank
(possibly a rational number, in case of ties) for each data item.

[list_begin arguments]
[arg_def list args] - Two or more lists of data
[list_end]
[para]

[call [cmd ::math::statistics::test-Wilcoxon] [arg sample_a] [arg sample_b]]
Compute the Wilcoxon test statistic to determine if two samples have the
same median or not. (The statistic can be regarded as standard normal, if the
sample sizes are both larger than 10. Returns the value of this statistic.

[list_begin arguments]
[arg_def list sample_a] - List of data comprising the first sample
[arg_def list sample_b] - List of data comprising the second sample
[list_end]
[para]

[call [cmd ::math::statistics::spearman-rank] [arg sample_a] [arg sample_b]]
Return the Spearman rank correlation as an alternative to the ordinary (Pearson's) correlation
coefficient. The two samples should have the same number of data.

[list_begin arguments]
[arg_def list sample_a] - First list of data
[arg_def list sample_b] - Second list of data
[list_end]
[para]

[call [cmd ::math::statistics::spearman-rank-extended] [arg sample_a] [arg sample_b]]
Return the Spearman rank correlation as an alternative to the ordinary (Pearson's) correlation
coefficient as well as additional data. The two samples should have the same number of data.
The procedure returns the correlation coefficient, the number of data pairs used and the
z-score, an approximately standard normal statistic, indicating the significance of the correlation.

[list_begin arguments]
[arg_def list sample_a] - First list of data
[arg_def list sample_b] - Second list of data
[list_end]

[call [cmd ::math::statistics::kernel-density] [arg data] opt [arg "-option value"] ...]]
Return the density function based on kernel density estimation. The procedure is controlled by
a small set of options, each of which is given a reasonable default.
[nl]
The return value consists of three lists: the centres of the bins, the associated probability
density and a list of computational parameters (begin and end of the interval, mean and standard
deviation and the used bandwidth). The computational parameters can be used for further analysis.

[list_begin arguments]
[arg_def list data] - The data to be examined
[arg_def list args] - Option-value pairs:
[list_begin definitions]
[def "[option -weights] [arg weights]"]  Per data point the weight (default: 1 for all data)
[def "[option -bandwidth] [arg value]"]  Bandwidth to be used for the estimation (default: determined from standard deviation)
[def "[option -number] [arg value]"]  Number of bins to be returned (default: 100)
[def "[option -interval] [arg "{begin end}"]"]  Begin and end of the interval for
which the density is returned (default: mean +/- 3*standard deviation)
[def "[option -kernel] [arg function]"]  Kernel to be used (One of: gaussian, cosine,
epanechnikov, uniform, triangular, biweight, logistic; default: gaussian)
[list_end]
[list_end]

[list_end]

[section "MULTIVARIATE LINEAR REGRESSION"]

Besides the linear regression with a single independent variable, the
statistics package provides two procedures for doing ordinary

[list_end]
[para]

[call [cmd ::math::statistics::subdivide]]
Routine [emph PM] - not implemented yet
[para]
































































[list_end]

[section "PLOT PROCEDURES"]
The following simple plotting procedures are available:
[list_begin definitions]

[call [cmd ::math::statistics::plot-scale] [arg canvas] \

#                Eric Kemp-Benedict, february 2007
# 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

package provide math::statistics 0.9
package require math

# ::math::statistics --
#   Namespace holding the procedures and variables
#

namespace eval ::math::statistics {

        if { $range < $rlower } { lappend result $i }
        if { $range > $rupper } { lappend result $i }
    }

    return $result
}





#
# 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]
source [file join [file dirname [info script]] mvlinreg.tcl]
source [file join [file dirname [info script]] kruskal.tcl]
source [file join [file dirname [info script]] wilcoxon.tcl]
source [file join [file dirname [info script]] stat_kernel.tcl]

#
# Define the tables
#
namespace eval ::math::statistics {
    variable student_t_table