[EKB] The ::math::statistics package in tcllib includes a linear model procedure for a single variable in which 

Y = a * X + b + error

However, sometimes there is a need to do multivariate regression of the form

Y = a1 * X1 + a2 * X2 + ... + aN * XN + b + error

This page has an implementation of multivariate linear regression. (I adapted it from my own [Substantively Weighted Least Squares (SWLS)] page.)

Here's the library (in the mvlinreg namespace):

======
    # Copyright 2007 Eric Kemp-Benedict
    # This library may be used and modified without notification to the author.
    # It would be greatly appreciated if corrections and improvements were
    # submitted to the appropriate page on the Tcler's Wiki.
    package require Tcl 8.4
    package require math::linearalgebra 1.0
    package require math::statistics 0.1.1
    
    # mvlinreg = Multivariate Linear Regression
    namespace eval mvlinreg {
        variable epsilon 1.0e-7
        
        namespace import -force \
            ::math::linearalgebra::mkMatrix \
            ::math::linearalgebra::mkVector \
            ::math::linearalgebra::mkIdentity \
            ::math::linearalgebra::mkDiagonal \
            ::math::linearalgebra::getrow \
            ::math::linearalgebra::setrow \
            ::math::linearalgebra::getcol \
            ::math::linearalgebra::setcol \
            ::math::linearalgebra::getelem \
            ::math::linearalgebra::setelem \
            ::math::linearalgebra::dotproduct \
            ::math::linearalgebra::matmul \
            ::math::linearalgebra::add \
            ::math::linearalgebra::sub \
            ::math::linearalgebra::solveGauss
        
        # NOTE: The authors of math::linearalgebra forgot to
        #   export ::math::linearalgebra::transpose
            
        ########################################
        #
        # t-stats
        #
        ########################################
        # mvlinreg::tstat n ?alpha?
        # Returns inverse of the single-tailed t distribution
        #  given number of degrees of freedom & confidence
        # Defaults to alpha = 0.05
        # Iterates until result is within epsilon of the target
        proc tstat {n {alpha 0.05}} {
            variable epsilon
            variable tvals
            
            if [info exists tvals($n:$alpha)] {
                return $tvals($n:$alpha)
            }
            
            set deltat [expr {100 * $epsilon}]
            # For one-tailed distribution, 
            set ptarg [expr {1.000 - $alpha/2.0}]
            set maxiter 100
            
            # Initial value for t
            set t 2.0
            
            set niter 0
            while {abs([::math::statistics::cdf-students-t $n $t] - $ptarg) > $epsilon} {
                set pstar [::math::statistics::cdf-students-t $n $t]
                set pl [::math::statistics::cdf-students-t $n [expr {$t - $deltat}]]
                set ph [::math::statistics::cdf-students-t $n [expr {$t + $deltat}]]
                
                set t [expr {$t + 2.0 * $deltat * ($ptarg - $pstar)/($ph - $pl)}]
                
                incr niter
                if {$niter == $maxiter} {
                    error "mvlinreg::tstat: Did not converge after $niter iterations"
                    return {}
                }
            }
            
            # Cache the result to shorten the call in future
            set tvals($n:$alpha) $t
            
            return $t
        }
            
        ########################################
        #
        # Weighted Least Squares
        #
        ########################################
        # mvlinreg::wls w [list y x's] w [list y x's] ...
        # Returns:
        #   R-squared
        #   Adj R-squared
        #   coefficients of x's in fit
        #   standard errors of the coefficients
        #   95% confidence bounds for coefficients
        proc wls {args} {
            
            # Fill the matrices of x & y values, and weights
            # For n points, k coefficients
            
            # The number of points is equal to half the arguments (n weights, n points)
            set n [expr {[llength $args]/2}]
            
            set firstloop true
            # Sum up all y values to take an average
            set ysum 0
            # Add up the weights
            set wtsum 0
            # Count over rows (points) as you go
            set point 0
            foreach {wt pt} $args {
                
                # Check inputs
                if {[string is double $wt] == 0} {
                    error "mvlinreg::wls: Weight \"$wt\" is not a number"
                    return {}
                }
                
                ## -- Check dimensions, initialize
                if $firstloop {
                    # k = num of vals in pt = 1 + number of x's (because of constant)
                    set k [llength $pt]
                    if {$n <= [expr {$k + 1}]} {
                        error "mvlinreg::wls: Too few degrees of freedom: $k variables but only $n points"
                        return {}
                    }
                    set X [mkMatrix $n $k]
                    set y [mkVector $n]
                    set I_x [mkIdentity $k]
                    set I_y [mkIdentity $n]
                    
                    set firstloop false
                    
                } else {
                    # Have to have same number of x's for all points
                    if {$k != [llength $pt]} {
                        error "mvlinreg::wls: Point \"$pt\" has wrong number of values (expected $k)"
                        # Clean up
                        return {}
                    }
                }
                
                ## -- Extract values from set of points
                # Make a list of y values
                set yval [expr {double([lindex $pt 0])}]
                setelem y $point $yval
                set ysum [expr {$ysum + $wt * $yval}]
                set wtsum [expr {$wtsum + $wt}]
                # Add x-values to the x-matrix
                set xrow [lrange $pt 1 end]
                # Add the constant (value = 1.0)
                lappend xrow 1.0
                setrow X $point $xrow
                # Create list of weights & square root of weights
                lappend w [expr {double($wt)}]
                lappend sqrtw [expr {sqrt(double($wt))}]
                
                incr point
                
            }
    
            set ymean [expr {double($ysum)/$wtsum}]
            set W [mkDiagonal $w]
            set sqrtW [mkDiagonal $sqrtw]
            
            # Calculate sum os square differences for x's
            for {set r 0} {$r < $k} {incr r} {
                set xsqrsum 0.0
                set xmeansum 0.0
                # Calculate sum of squared differences as: sum(x^2) - (sum x)^2/n
                for {set t 0} {$t < $n} {incr t} {
                    set xval [getelem $X $t $r]
                    set xmeansum [expr {$xmeansum + double($xval)}]
                    set xsqrsum [expr {$xsqrsum + double($xval * $xval)}]
                }
                lappend xsqr [expr {$xsqrsum - $xmeansum * $xmeansum/$n}]
            }
            
            ## -- Set up the X'WX matrix
            set XtW [matmul [::math::linearalgebra::transpose $X] $W]
            set XtWX [matmul $XtW $X]
            
            # Invert
            set M [solveGauss $XtWX $I_x]
            
            set beta [matmul $M [matmul $XtW $y]]
            
            ### -- Residuals & R-squared
            # 1) Generate list of diagonals of the hat matrix
            set H [matmul $X [matmul $M $XtW]]
            for {set i 0} {$i < $n} {incr i} {
                lappend h_ii [getelem $H $i $i]
            }
    
            set R [matmul $sqrtW [matmul [sub $I_y $H] $y]]
            set yhat [matmul $H $y]
            
            # 2) Generate list of residuals, sum of squared residuals, r-squared
            set sstot 0.0
            set ssreg 0.0
            # Note: Relying on representation of Vector as a list for y, yhat
            foreach yval $y wt $w yhatval $yhat {
                set sstot [expr {$sstot + $wt * ($yval - $ymean) * ($yval - $ymean)}]
                set ssreg [expr {$ssreg + $wt * ($yhatval - $ymean) * ($yhatval - $ymean)}]
            }
            set r2 [expr {double($ssreg)/$sstot}]
            set adjr2 [expr {1.0 - (1.0 - $r2) * ($n - 1)/($n - $k)}]
            set sumsqresid [dotproduct $R $R]
            set s2 [expr {double($sumsqresid) / double($n - $k)}]
            
            ### -- Confidence intervals for coefficients
            set tvalue [tstat [expr {$n - $k}]]
            for {set i 0} {$i < $k} {incr i} {
                set stderr [expr {sqrt($s2 * [getelem $M $i $i])}]
                set mid [lindex $beta $i]
                lappend stderrs $stderr
                lappend confinterval [list [expr {$mid - $tvalue * $stderr}] [expr {$mid + $tvalue * $stderr}]]
            }
    
            return [list $r2 $adjr2 $beta $stderrs $confinterval]
        }
        
        ########################################
        #
        # Ordinary (unweighted) Least Squares
        #
        ########################################
        # mvlinreg::ols [list y x's] [list y x's] ...
        # Returns:
        #   R-squared
        #   Adj R-squared
        #   coefficients of x's in fit
        #   standard errors of the coefficients
        #   95% confidence bounds for coefficients
        proc ols {args} {
            set newargs {}
            foreach pt $args {
                lappend newargs 1 $pt
            }
            return [eval wls $newargs]
        }
    }

if 0 {
Here's an example:
}

    
    ##
    ##
    ## TEST IT
    ##
    ##
    
    set data {{183.34 100. 9.43} \
    {222.67 100. 9.14} \
    {410.77 20.3 8.81} \
    {268.57 17.53 8.74} \
    {499.28 10.95 8.71} \
    {385.97 51.81 8.78} \
    {474.53 17.01 8.74} \
    {550.04 16.6 8.72} \
    {221.54 100. 9.49} \
    {284.93 49.9 9.01} \
    {324.19 23.02 9.37} \
    {364.6 16.59 8.7} \
    {338.62 20.03 9.13} \
    {543.19 13.44 8.75} \
    {356.96 25.64 8.94} \
    {296.54 18.63 8.74} \
    {336.98 50.07 8.82} \
    {398.5 15.35 8.85} \
    {216.37 47.41 9.32} \
    {325.02 14.83 8.9} \
    {315.55 9.43 8.83} \
    {321.57 76.98 8.89} \
    {443.21 51.28 8.8} \
    {493.22 13.19 8.69} \
    {554.79 19.62 8.89} \
    {527.46 63.23 8.75} \
    {490.59 10.72 8.72} \
    {389.31 14.39 8.73} \
    {218.9 8.65 8.71} \
    {315.65 7.29 8.84}}
    
    set results [eval mvlinreg::ols $data]
    
    puts "R-squared: [lindex $results 0]"
    puts "Adj R-squared: [lindex $results 1]"
    puts "Coefficients \u00B1 s.e. -- \[95% confidence interval\]:"
    foreach val [lindex $results 2] se [lindex $results 3] bounds [lindex $results 4] {
        set lb [lindex $bounds 0]
        set ub [lindex $bounds 1]
        puts "   $val \u00B1 $se -- \[$lb to $ub\]"
    }
    
if 0 {
This should give the following output:
 R-squared: 0.371918955025
 Adj R-squared: 0.325394433175
 Coefficients ± s.e. -- [95% confidence interval]:
   -0.432215837601 ± 0.744774375289 -- [-1.96036662835 to 1.09593495315]
   -250.868151309 ± 92.4723901677 -- [-440.605823342 to -61.1304792767]
   2615.78338226 ± 807.559562552 -- [958.808028331 to 4272.75873618]
}
======


----
[AM] (19 february 2007) You are welcome to contribute this to Tcllib :) Just say the word, and I will incorporate 
it in the package. 

[EKB] I'd be delighted to have it in Tcllib! What should I do to make that happen?

[LV] One place developers can start is visiting http://tcllib.sf.net/
and posting a new Feature Request, specifying that you are the creator
of the code in question, specifying the license you use (tcllib uses
BSD based licensing, I believe), and mention that we'd asked you 
to contribute it. Man pages, test cases, demos, etc. are all wonderful
things to have available, if possible. But if they are not available
and you don't have time/interest in generating them, then perhaps someone
who is a fan of the module could write them.

[EKB] Thanks. Sounds good. I'll work on creating demos and a man page.

[EKB] Done! I've submitted it, with a plain text man page (no *roff) to tcllib at sf.net. I hope I've included what's needed!

<<categories>> Mathematics | Statistics