[manpage_begin math::statistics n 0.2]
[moddesc {Tcl Math Library}]
[titledesc {Basic statistical functions and procedures}]
[require Tcl 8]
[require math::statistics 0.2]
[description]
[para]

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

[list_begin bullet]
[bullet]
Descriptive statistical parameters (mean, minimum, maximum, standard
deviation)

[bullet]
Estimates of the distribution in the form of histograms and quantiles

[bullet]
Basic testing of hypotheses

[bullet]
Probability and cumulative density functions

[list_end]
It is meant to help in developing data analysis applications or doing
ad hoc data analysis, it is not in itself a full application, nor is it
intended to rival with full (non-)commercial statistical packages.

[para]
The purpose of this document is to describe the implemented procedures
and provide some examples of their usage. As there is ample literature
on the algorithms involved, we refer to relevant text books for more
explanations.

The package contains a fairly large number of public procedures. They
can be distinguished in three sets: general procedures, procedures
that deal with specific statistical distributions, list procedures to
select or transform data and simple plotting procedures (these require
Tk).

[emph Note:] The data that need to be analyzed are always contained in a
simple list. Missing values are represented as empty list elements.

[section "GENERAL PROCEDURES"]
The general statistical procedures are:
[list_begin definitions]

[call [cmd ::math::statistics::mean] [arg data]]
Determine the [emph mean] value of the given list of data.
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::min] [arg data]]
Determine the [emph minimum] value of the given list of data.
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::max] [arg data]]
Determine the [emph maximum] value of the given list of data.
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::number] [arg data]]
Determine the [emph number] of non-missing data in the given list
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::stdev] [arg data]]
Determine the [emph "standard deviation"] of the data in the given list
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::var] [arg data]]
Determine the [emph variance] of the data in the given list
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::median] [arg data]]
Determine the [emph median] of the data in the given list
(Note that this requires sorting the data, which may be a
costly operation)
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::basic-stats] [arg data]]
Determine a list of all the descriptive parameters: mean, minimum,
maximum, number of data, standard deviation and variance.
[nl]
(This routine is called whenever either or all of the basic statistical
parameters are required. Hence all calculations are done and the
relevant values are returned.)
[nl]
[arg data] - List of data
[nl]

[call [cmd ::math::statistics::histogram] [arg limits] [arg values]]
Determine histogram information for the given list of data. Returns a
list consisting of the number of values that fall into each interval.
(The first interval consists of all values lower than the first limit,
the last interval consists of all values greater than the last limit.
There is one more interval than there are limits.)
[nl]
[arg limits] - List of upper limits (in ascending order) for the
intervals of the histogram.
[nl]
[arg values] - List of data
[nl]

[call [cmd ::math::statistics::corr] [arg data1] [arg data2]]
Determine the correlation coefficient between two sets of data.
[nl]
[arg data1] - First list of data
[nl]
[arg data2] - Second list of data
[nl]

[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)
[nl]
[arg data] - List of raw data values (small sample)
[nl]
[arg confidence] - Confidence level (0.95 or 0.99 for instance)
[nl]

[call [cmd ::math::statistics::t-test-mean] [arg data] [arg est_mean] \
[arg est_stdev] [arg confidence]]
Test whether the mean value of a sample is in accordance with the
estimated normal distribution with a certain level of confidence.
Returns 1 if the test succeeds or 0 if the mean is unlikely to fit
the given distribution.
[nl]
[arg data] - List of raw data values (small sample)
[nl]
[arg est_mean] - Estimated mean of the distribution
[nl]
[arg est_stdev] - Estimated stdev of the distribution
[nl]
[arg confidence] - Confidence level (0.95 or 0.99 for instance)
[nl]

[call [cmd ::math::statistics::quantiles] [arg data] [arg confidence]]
Return the quantiles for a given set of data
[nl]
[arg data] - List of raw data values
[nl]
[arg confidence] - Confidence level (0.95 or 0.99 for instance)
[nl]

[call [cmd ::math::statistics::quantiles] [arg limits] [arg counts] [arg confidence]]
Return the quantiles based on histogram information (alternative to the
call with two arguments)
[nl]
[arg limits] - List of upper limits from histogram
[nl]
[arg counts] - List of counts for for each interval in histogram
[nl]
[arg confidence] -  Confidence level (0.95 or 0.99 for instance)
[nl]

[call [cmd ::math::statistics::autocorr] [arg data]]
Return the autocorrelation function as a list of values (assuming
equidistance between samples, about 1/2 of the number of raw data)
[nl]
The correlation is determined in such a way that the first value is
always 1 and all others are equal to or smaller than 1. The number of
values involved will diminish as the "time" (the index in the list of
returned values) increases
[nl]
[arg data] - Raw data for which the autocorrelation must be determined
[nl]

[call [cmd ::math::statistics::crosscorr] [arg data1] [arg data2]]
Return the cross-correlation function as a list of values (assuming
equidistance between samples, about 1/2 of the number of raw data)
[nl]
The correlation is determined in such a way that the values can never
exceed 1 in magnitude. The number of values involved will diminish
as the "time" (the index in the list of returned values) increases.
[nl]
[arg data1] - First list of data
[nl]
[arg data2] - Second list of data
[nl]

[call [cmd ::math::statistics::mean-histogram-limits] [arg mean] \
[arg stdev] [arg number]]
Determine reasonable limits based on mean and standard deviation
for a histogram
[nl]
Convenience function - the result is suitable for the histogram function.
[nl]
[arg mean] - Mean of the data
[nl]
[arg stdev] - Standard deviation
[nl]
[arg number] - Number of limits to generate (defaults to 8)
[nl]

[call [cmd ::math::statistics::minmax-histogram-limits] [arg min] \
[arg max] [arg number]]
Determine reasonable limits based on a minimum and maximum for a histogram
[nl]
Convenience function - the result is suitable for the histogram function.
[nl]
[arg min] - Expected minimum
[nl]
[arg max] - Expected maximum
[nl]
[arg number] - Number of limits to generate (defaults to 8)
[nl]

[call [cmd ::math::statistics::linear-model] [arg xdata] \
[arg ydata] [arg intercept]]
Determine the coefficients for a linear regression between
two series of data (the model: Y = A + B*X). Returns a list of
parameters describing the fit
[nl]
[arg xdata] - List of independent data
[nl]
[arg ydata] - List of dependent data to be fitted
[nl]
[arg intercept] - (Optional) compute the intercept (1, default) or fit
to a line through the origin (0)
[nl]
The result consists of the following list:
[list_begin bullet]
[bullet]
(Estimate of) Intercept A
[bullet]
(Estimate of) Slope B
[bullet]
Standard deviation of Y relative to fit
[bullet]
Correlation coefficient R2
[bullet]
Number of degrees of freedom df
[bullet]
Standard error of the intercept A
[bullet]
Significance level of A
[bullet]
Standard error of the slope B
[bullet]
Significance level of B
[list_end]

[call [cmd ::math::statistics::linear-residuals] [arg xdata] [arg ydata] \
[arg intercept]]
Determine the difference between actual data and predicted from
the linear model.
[nl]
Returns a list of the differences between the actual data and the
predicted values.
[nl]
[arg xdata] - List of independent data
[nl]
[arg ydata] - List of dependent data to be fitted
[nl]
[arg intercept] - (Optional) compute the intercept (1, default) or fit
to a line through the origin (0)

[call [cmd ::math::statistics::test-2x2] [arg n11] [arg n21] [arg n12] [arg n22]]
Determine if two set of samples, each from a binomial distribution,
differ significantly or not (implying a different parameter).
[nl]
Returns the "chi-square" value, which can be used to the determine the
significance.
[nl]
[arg n11] - Number of outcomes with the first value from the first sample.
[nl]
[arg n21] - Number of outcomes with the first value from the second sample.
[nl]
[arg n12] - Number of outcomes with the second value from the first sample.
[nl]
[arg n22] - Number of outcomes with the second value from the second sample.
[nl]

[call [cmd ::math::statistics::print-2x2] [arg n11] [arg n21] [arg n12] [arg n22]]
Determine if two set of samples, each from a binomial distribution,
differ significantly or not (implying a different parameter).
[nl]
Returns a short report, useful in an interactive session.
[nl]
[arg n11] - Number of outcomes with the first value from the first sample.
[nl]
[arg n21] - Number of outcomes with the first value from the second sample.
[nl]
[arg n12] - Number of outcomes with the second value from the first sample.
[nl]
[arg n22] - Number of outcomes with the second value from the second sample.
[nl]

[call [cmd ::math::statistics::control-xbar] [arg data] [opt nsamples]]
Determine the control limits for an xbar chart. The number of data
in each subsample defaults to 4. At least 20 subsamples are required.
[nl]
Returns the mean, the lower limit, the upper limit and the number of
data per subsample.
[nl]
[arg data] - List of observed data
[nl]
[arg nsamples] - Number of data per subsample
[nl]

[call [cmd ::math::statistics::control-Rchart] [arg data] [opt nsamples]]
Determine the control limits for an R chart. The number of data
in each subsample defaults to 4. At least 20 subsamples are required.
[nl]
Returns the mean range, the lower limit, the upper limit and the number
of data per subsample.
[nl]
[arg data] - List of observed data
[nl]
[arg nsamples] - Number of data per subsample
[nl]

[call [cmd ::math::statistics::test-xbar] [arg control] [arg data]]
Determine if the data exceed the control limits for the xbar chart.
[nl]
Returns a list of subsamples (their indices) that indeed violate the
limits.
[nl]
[arg control] - Control limits as returned by the "control-xbar" procedure
[nl]
[arg data] - List of observed data
[nl]

[call [cmd ::math::statistics::test-Rchart] [arg control] [arg data]]
Determine if the data exceed the control limits for the R chart.
[nl]
Returns a list of subsamples (their indices) that indeed violate the
limits.
[nl]
[arg control] - Control limits as returned by the "control-Rchart" procedure
[nl]
[arg data] - List of observed data
[nl]

[list_end]

[section "STATISTICAL DISTRIBUTIONS"]
In the literature a large number of probability distributions can be
found. The statistics package supports:
[list_begin bullet]
[bullet]
The normal or Gaussian distribution
[bullet]
The uniform distribution - equal probability for all data within a given
interval
[bullet]
The exponential distribution - useful as a model for certain
extreme-value distributions.
[bullet]
PM - binomial, Poisson, chi-squared, student's T, F.
[list_end]

In principle for each distribution one has procedures for:
[list_begin bullet]
[bullet]
The probability density (pdf-*)
[bullet]
The cumulative density (cdf-*)
[bullet]
Quantiles for the given distribution (quantiles-*)
[bullet]
Histograms for the given distribution (histogram-*)
[bullet]
List of random values with the given distribution (random-*)
[list_end]

The following procedures have been implemented:
[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.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg stdev] - Standard deviation of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[call [cmd ::math::statistics::pdf-exponential] [arg mean] [arg value]]
Return the probability of a given value for an exponential
distribution with given mean.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[call [cmd ::math::statistics::pdf-uniform] [arg xmin] [arg xmax] [arg value]]
Return the probability of a given value for a uniform
distribution with given extremes.
[nl]
[arg xmin] - Minimum value of the distribution
[nl]
[arg xmin] - Maximum value of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[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.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg stdev] - Standard deviation of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[call [cmd ::math::statistics::cdf-exponential] [arg mean] [arg value]]
Return the cumulative probability of a given value for an exponential
distribution with given mean.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[call [cmd ::math::statistics::cdf-uniform] [arg xmin] [arg xmax] [arg value]]
Return the cumulative probability of a given value for a uniform
distribution with given extremes.
[nl]
[arg xmin] - Minimum value of the distribution
[nl]
[arg xmin] - Maximum value of the distribution
[nl]
[arg value] - Value for which the probability is required
[nl]

[call [cmd ::math::statistics::cdf-students-t] [arg degrees] [arg value]]
Return the cumulative probability of a given value for a Student's t
distribution with given number of degrees.
[nl]
[arg degrees] - Number of degrees of freedom
[nl]
[arg value] - Value for which the probability is required
[nl]

[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.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg stdev] - Standard deviation of the distribution
[nl]
[arg number] - Number of values to be returned
[nl]

[call [cmd ::math::statistics::random-exponential] [arg mean] [arg number]]
Return a list of "number" random values satisfying an exponential
distribution with given mean.
[nl]
[arg mean] - Mean value of the distribution
[nl]
[arg number] - Number of values to be returned
[nl]

[call [cmd ::math::statistics::random-uniform] [arg xmin] [arg xmax] [arg value]]
Return a list of "number" random values satisfying a uniform
distribution with given extremes.
[nl]
[arg xmin] - Minimum value of the distribution
[nl]
[arg xmin] - Maximum value of the distribution
[nl]
[arg number] - Number of values to be returned
[nl]

[call [cmd ::math::statistics::histogram-uniform] [arg xmin] [arg xmax] [arg limits] [arg number]]
Return the expected histogram for a uniform distribution.
[nl]
[arg xmin] - Minimum value of the distribution
[nl]
[arg xmax] - Maximum value of the distribution
[nl]
[arg limits] - Upper limits for the buckets in the histogram
[nl]
[arg number] - Total number of "observations" in the histogram
[nl]


[list_end]
TO DO: more function descriptions to be added

[section "DATA MANIPULATION"]
The data manipulation procedures act on lists or lists of lists:

[list_begin definitions]

[call [cmd ::math::statistics::filter] [arg varname] [arg data] [arg expression]]
Return a list consisting of the data for which the logical
expression is true (this command works analogously to the command [cmd foreach]).
[nl]
[arg varname] - Name of the variable used in the expression
[nl]
[arg data] - List of data
[nl]
[arg expression] - Logical expression using the variable name
[nl]

[call [cmd ::math::statistics::map] [arg varname] [arg data] [arg expression]]
Return a list consisting of the data that are transformed via the
expression.
[nl]
[arg varname] - Name of the variable used in the expression
[nl]
[arg data] - List of data
[nl]
[arg expression] - Expression to be used to transform (map) the data
[nl]

[call [cmd ::math::statistics::samplescount] [arg varname] [arg list] [arg expression]]
Return a list consisting of the [emph counts] of all data in the
sublists of the "list" argument for which the expression is true.
[nl]
[arg varname] - Name of the variable used in the expression
[nl]
[arg data] - List of sublists, each containing the data
[nl]
[arg expression] - Logical expression to test the data (defaults to
"true").
[nl]

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

[list_end]

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

[call [cmd ::math::statistics::plot-scale] [arg canvas] \
[arg xmin] [arg xmax] [arg ymin] [arg ymax]]
Set the scale for a plot in the given canvas. All plot routines expect
this function to be called first. There is no automatic scaling
provided.
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg xmin] - Minimum x value
[nl]
[arg xmax] - Maximum x value
[nl]
[arg ymin] - Minimum y value
[nl]
[arg ymax] - Maximum y value
[nl]

[call [cmd ::math::statistics::plot-xydata] [arg canvas] \
[arg xdata] [arg ydata] [arg tag]]
Create a simple XY plot in the given canvas - the data are
shown as a collection of dots. The tag can be used to manipulate the
appearance.
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg xdata] - Series of independent data
[nl]
[arg ydata] - Series of dependent data
[nl]
[arg tag] - Tag to give to the plotted data (defaults to xyplot)
[nl]

[call [cmd ::math::statistics::plot-xyline] [arg canvas] \
[arg xdata] [arg ydata] [arg tag]]
Create a simple XY plot in the given canvas - the data are
shown as a line through the data points. The tag can be used to
manipulate the appearance.
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg xdata] - Series of independent data
[nl]
[arg ydata] - Series of dependent data
[nl]
[arg tag] - Tag to give to the plotted data (defaults to xyplot)
[nl]

[call [cmd ::math::statistics::plot-tdata] [arg canvas] \
[arg tdata] [arg tag]]
Create a simple XY plot in the given canvas - the data are
shown as a collection of dots. The horizontal coordinate is equal to the
index. The tag can be used to manipulate the appearance.
This type of presentation is suitable for autocorrelation functions for
instance or for inspecting the time-dependent behaviour.
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg tdata] - Series of dependent data
[nl]
[arg tag] - Tag to give to the plotted data (defaults to xyplot)
[nl]

[call [cmd ::math::statistics::plot-tline] [arg canvas] \
[arg tdata] [arg tag]]
Create a simple XY plot in the given canvas - the data are
shown as a line. See plot-tdata for an explanation.
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg tdata] - Series of dependent data
[nl]
[arg tag] - Tag to give to the plotted data (defaults to xyplot)
[nl]

[call [cmd ::math::statistics::plot-histogram] [arg canvas] \
[arg counts] [arg limits] [arg tag]]
Create a simple histogram in the given canvas
[nl]
[arg canvas] - Canvas widget to use
[nl]
[arg counts] - Series of bucket counts
[nl]
[arg limits] - Series of upper limits for the buckets
[nl]
[arg tag] - Tag to give to the plotted data (defaults to xyplot)
[nl]

[list_end]

[section {THINGS TO DO}]
The following procedures are yet to be implemented:
[list_begin bullet]
[bullet]
F-test-stdev
[bullet]
interval-mean-stdev
[bullet]
histogram-normal
[bullet]
histogram-exponential
[bullet]
test-histogram
[bullet]
test-corr
[bullet]
quantiles-*
[bullet]
fourier-coeffs
[bullet]
fourier-residuals
[bullet]
onepar-function-fit
[bullet]
onepar-function-residuals
[bullet]
plot-linear-model
[bullet]
subdivide
[list_end]

[section EXAMPLES]
The code below is a small example of how you can examine a set of
data:
[para]
[example_begin]

# Simple example:
# - Generate data (as a cheap way of getting some)
# - Perform statistical analysis to describe the data
#
package require math::statistics

#
# Two auxiliary procs
#
proc pause {time} {
   set wait 0
   after [lb]expr {$time*1000}[rb] {set ::wait 1}
   vwait wait
}

proc print-histogram {counts limits} {
   foreach count $counts limit $limits {
      if { $limit != {} } {
         puts [lb]format "<%12.4g\t%d" $limit $count[rb]
         set prev_limit $limit
      } else {
         puts [lb]format ">%12.4g\t%d" $prev_limit $count[rb]
      }
   }
}

#
# Our source of arbitrary data
#
proc generateData { data1 data2 } {
   upvar 1 $data1 _data1
   upvar 1 $data2 _data2

   set d1 0.0
   set d2 0.0
   for { set i 0 } { $i < 100 } { incr i } {
      set d1 [lb]expr {10.0-2.0*cos(2.0*3.1415926*$i/24.0)+3.5*rand()}[rb]
      set d2 [lb]expr {0.7*$d2+0.3*$d1+0.7*rand()}[rb]
      lappend _data1 $d1
      lappend _data2 $d2
   }
   return {}
}

#
# The analysis session
#
package require Tk
console show
canvas .plot1
canvas .plot2
pack   .plot1 .plot2 -fill both -side top

generateData data1 data2

puts "Basic statistics:"
set b1 [lb]::math::statistics::basic-stats $data1[rb]
set b2 [lb]::math::statistics::basic-stats $data2[rb]
foreach label {mean min max number stdev var} v1 $b1 v2 $b2 {
   puts "$label\t$v1\t$v2"
}
puts "Plot the data as function of \"time\" and against each other"
::math::statistics::plot-scale .plot1  0 100  0 20
::math::statistics::plot-scale .plot2  0 20   0 20
::math::statistics::plot-tline .plot1 $data1
::math::statistics::plot-tline .plot1 $data2
::math::statistics::plot-xydata .plot2 $data1 $data2

puts "Correlation coefficient:"
puts [lb]::math::statistics::corr $data1 $data2]

pause 2
puts "Plot histograms"
.plot2 delete all
::math::statistics::plot-scale .plot2  0 20 0 100
set limits         [lb]::math::statistics::minmax-histogram-limits 7 16[rb]
set histogram_data [lb]::math::statistics::histogram $limits $data1[rb]
::math::statistics::plot-histogram .plot2 $histogram_data $limits

puts "First series:"
print-histogram $histogram_data $limits

pause 2
set limits         [lb]::math::statistics::minmax-histogram-limits 0 15 10[rb]
set histogram_data [lb]::math::statistics::histogram $limits $data2[rb]
::math::statistics::plot-histogram .plot2 $histogram_data $limits d2
.plot2 itemconfigure d2 -fill red

puts "Second series:"
print-histogram $histogram_data $limits

puts "Autocorrelation function:"
set  autoc [lb]::math::statistics::autocorr $data1[rb]
puts [lb]::math::statistics::map $autoc {[lb]format "%.2f" $x]}[rb]
puts "Cross-correlation function:"
set  crossc [lb]::math::statistics::crosscorr $data1 $data2[rb]
puts [lb]::math::statistics::map $crossc {[lb]format "%.2f" $x[rb]}[rb]

::math::statistics::plot-scale .plot1  0 100 -1  4
::math::statistics::plot-tline .plot1  $autoc "autoc"
::math::statistics::plot-tline .plot1  $crossc "crossc"
.plot1 itemconfigure autoc  -fill green
.plot1 itemconfigure crossc -fill yellow

puts "Quantiles: 0.1, 0.2, 0.5, 0.8, 0.9"
puts "First:  [lb]::math::statistics::quantiles $data1 {0.1 0.2 0.5 0.8 0.9}[rb]"
puts "Second: [lb]::math::statistics::quantiles $data2 {0.1 0.2 0.5 0.8 0.9}[rb]"

[example_end]
If you run this example, then the following should be clear:
[list_begin bullet]
[bullet]
There is a strong correlation between two time series, as displayed by
the raw data and especially by the correlation functions.
[bullet]
Both time series show a significant periodic component
[bullet]
The histograms are not very useful in identifying the nature of the time
series - they do not show the periodic nature.
[list_end]

[keywords mathematics "data analysis" statistics]

[manpage_end]