# pca_example.tcl --
#     Follow the PCA example from NIST (as paraphrased by Ed Hume, using the
#     LA package), employing the math::linearalgebra package
#
#     Note:
#     I may add a few procs to the linearalgebra package
#
package require math::linearalgebra

# Sample data - from the NIST page on the subject
set data {
    {7 4 3}
    {4 1 8}
    {6 3 5}
    {8 6 1}
    {8 5 7}
    {7 2 9}
    {5 3 3}
    {9 5 8}
    {7 4 5}
    {8 2 2}
}

# Normalise the matrix

set z [::math::linearalgebra::normalizeStat $data]

# Perform a singular value decomposition

set usv [::math::linearalgebra::determineSVD $z]

# Print the singular values

set eigenVectors [lindex $usv 2]
set singular     [lindex $usv 1]
puts "Singular values: $singular"

# Square the vector of singular values and divide by the number of observations - 1
# to get the eigenvalues

set factor [expr {1.0/([llength $data] - 1)}]

set eigen {}
foreach c $singular {
    lappend eigen [expr {$c**2 * $factor}]
}

puts "Eigenvalues: $eigen"
puts "Eigenvectors:\n[::math::linearalgebra::show $eigenVectors]"

puts "Verify that V*Vt = I"
puts "[::math::linearalgebra::show [::math::linearalgebra::matmul $eigenVectors [::math::linearalgebra::transpose $eigenVectors]]]"

# Determine the square roots of the eigenvalues
# We need them for the factor structure matrix
# At the same time, construct the diagonal with the reciprocals
set sqeigen {}
set rsqeigen {}
foreach c $eigen {
    lappend sqeigen  [expr {sqrt($c)}]
    lappend rsqeigen [expr {1.0/sqrt($c)}]
}

puts "Square root: $sqeigen"
set sqdiagonal [::math::linearalgebra::mkDiagonal $sqeigen]

# Now the factor structure matrix
set s [::math::linearalgebra::matmul $eigenVectors $sqdiagonal]

puts "Factor structure:\n[::math::linearalgebra::show $s]"

# Retain the first two factors only and print the communality
set s2 {}
foreach row $s {
    lappend s2 [lrange $row 0 1]
}

set s2s2 [::math::linearalgebra::matmul $s2 [::math::linearalgebra::transpose $s2]]
puts "SSt - first two factors\n[::math::linearalgebra::show $s2s2]"

set communality {}
for {set i 0} {$i < [llength $s2s2]} {incr i} {
    lappend communality [lindex $s2s2 $i $i]
}
puts "Communality: $communality"

# Now construct the factor score coefficients matrix (b)
set b [::math::linearalgebra::matmul $eigenVectors [::math::linearalgebra::mkDiagonal $rsqeigen]]

puts "Factor score coefficients:\n[::math::linearalgebra::show $b]"

# Examine the first observation - factor scores via the eigenvectors
set z1 [lindex $z 0]

set t [::math::linearalgebra::matmul $z1 $eigenVectors]
puts "Scores: $t"

set t2 0.0
foreach c $t e $eigen {
    set t2 [expr {$t2 + $c**2 / $e}]
}
puts "Hotelling: $t2"

# Fit the first observation, z1, with two principal components
set v2 {}
foreach row $eigenVectors {
    lappend v2 [lrange $row 0 1]
}
set t    [::math::linearalgebra::matmul $z1 $v2]
set zhat [::math::linearalgebra::matmul $t [::math::linearalgebra::transpose $v2]]

puts "z1:\n[::math::linearalgebra::show $z1]"
puts "zhat:\n[::math::linearalgebra::show $zhat]"
set difference  [::math::linearalgebra::show [::math::linearalgebra::sub $z1 $zhat]]
puts "difference:\n[::math::linearalgebra::show $difference]"
set qstat [::math::linearalgebra::dotproduct $difference $difference]
puts "Q statistic: $qstat"
                                                                   # Number of components
set qcorr [expr {$qstat * [llength $data] / double([llength $data] - 2 - 1)}]
puts "Q corrected: $qcorr"