stat: add matrix conditioning functions for MCA and MDS

WIP

stat/pca_helpers.go

@@ -0,0 +1,101 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stat
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/mat"
+)
+
+// TransformCDT transforms a Complete Disjunctive Table to a Transformed
+// Complete Disjunctive Table in order to be able to perform Multiple
+// Correspondence Analysis using PC.PrincipalComponents. TransformCDT places
+// the TCDT in dst and returns it.
+//
+// If dst is nil, a new mat.Dense is allocated. If dst is not a zero matrix,
+// the dimensions of dst and cdt must match otherwise TransformCDT will panic.
+// If cdt contains values other 0 or 1 TransformCDT will panic.
+//
+// It is safe to reuse cdt as dst.
+func TransformCDT(dst *mat.Dense, cdt mat.Matrix) *mat.Dense {
+	r, c := cdt.Dims()
+
+	if dst == nil {
+		dst = mat.NewDense(r, c, nil)
+	} else if dr, dc := dst.Dims(); !dst.IsZero() && (dr != r || dc != c) {
+		panic(mat.ErrShape)
+	}
+
+	for j := 0; j < c; j++ {
+		var p float64
+		for i := 0; i < r; i++ {
+			v := cdt.At(i, j)
+			if v != 0 && v != 1 {
+				panic("stat: input is not a complete disjunctive table")
+			}
+			p += v
+		}
+		for i := 0; i < r; i++ {
+			dst.Set(i, j, cdt.At(i, j)/p-1)
+		}
+	}
+
+	return dst
+}
+
+// ClassicalScaling converts a dissimilarity matrix to a matrix containing
+// Euclidean coordinates in order to be able to perform Torgerson's Classical
+// Multidimensional Scaling using PC.PrincipleComponenets. ClassicalScaling places
+// the coordinates in dst and returns it.
+//
+// If dst is nil, a new mat.Dense is allocated. If dst is not a zero matrix,
+// the dimensions of dst and dis must match otherwise ClassicalScaling will panic.
+// The dis matrix must be square or ClassicalScaling will panic.
+//
+// It is safe to reuse dis as dst.
+func ClassicalScaling(dst *mat.Dense, dis mat.Matrix) *mat.Dense {
+	// https://doi.org/10.1007/0-387-28981-X_12
+
+	r, c := dis.Dims()
+	if r != c {
+		panic(mat.ErrShape)
+	}
+
+	if dst == nil {
+		dst = mat.NewDense(r, c, nil)
+	} else if dr, dc := dst.Dims(); !dst.IsZero() && (dr != r || dc != c) {
+		panic(mat.ErrShape)
+	}
+
+	var gmean float64
+	rmean := make([]float64, r)
+	cmean := make([]float64, c)
+	dst.Apply(func(i, j int, v float64) float64 {
+		v *= v
+		gmean += v
+		rmean[i] += v
+		cmean[c] += v
+		return v
+	}, dis)
+	gmean /= float64(r * c)
+	floats.Scale(float64(c), rmean)
+	floats.Scale(float64(r), cmean)
+	dst.Apply(func(i, j int, v float64) float64 {
+		return -0.5 * (v - (rmean[i] + cmean[j]) + gmean)
+	}, dst)
+
+	var svd mat.SVD
+	svd.Factorize(dst, mat.SVDThin)
+	svd.UTo(dst)
+	vec := svd.Values(nil)
+	for i, v := range vec {
+		vec[i] = math.Sqrt(v)
+	}
+	dst.Mul(dst, mat.NewVecDense(len(vec), vec))
+
+	return dst
+}