# 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"

Category Mathematics | Category Statistics | [Category Data Analysis] |