Add CorrelationMatrix and Cov2Corr functions

This change includes a function to calculate correlation the
correlation matrix of input data, and a function which can convert
covariance matrices to correlation matrices.

I will be adding more test cases.

covariancematrix.go

@@ -80,3 +80,43 @@ func CovarianceMatrix(cov *mat64.Dense, x mat64.Matrix, wts []float64) *mat64.De
 	cov.Scale(1/(n-1), cov)
 	return cov
 }
+
+// CorrelationMatrix calculates a correlation matrix from a matrix of data,
+// using a two-pass algorithm. The matrix returned will be symmetric, square,
+// and positive-semidefinite.
+//
+// The weights wts should have the length equal to the number of rows in
+// input data matrix x. cov should either be a square matrix with the same
+// number of columns as the input data matrix x, or nil in which case a new
+// Dense matrix will be constructed.  Weights cannot be negative.
+func CorrelationMatrix(corr *mat64.Dense, x mat64.Matrix, wts []float64) *mat64.Dense {
+
+	// TODO(jonlawlor): indicate that the resulting matrix is symmetric, and change
+	// the returned type from a *mat.Dense to a *mat.Symmetric.
+
+	// This will panic if the sizes don't match, or if wts is the wrong size.
+	corr = CovarianceMatrix(corr, x, wts)
+	Cov2Corr(corr)
+	return corr
+}
+
+// Cov2Corr converts a covariance matrix to a correlation matrix.
+func Cov2Corr(cov *mat64.Dense) {
+
+	// TODO(jonlawlor): use a *mat64.Symmetric as input.
+
+	r, _ := cov.Dims()
+
+	s := make([]float64, r*r)
+	for i := 0; i < r; i++ {
+		s[i*(r+1)] = 1 / math.Sqrt(cov.At(i, i))
+	}
+	sMat := mat64.NewDense(r, r, s)
+	cov.Mul(sMat, cov)
+	cov.Mul(cov, sMat)
+
+	// Ensure that the diagonal is exactly ones.
+	for i := 0; i < r; i++ {
+		cov.Set(i, i, 1)
+	}
+}

covariancematrix_test.go

@@ -96,6 +96,71 @@ func TestCovarianceMatrix(t *testing.T) {
 	}
 }
 
+func TestCorrelationMatrix(t *testing.T) {
+	for i, test := range []struct {
+		data    *mat64.Dense
+		weights []float64
+		ans     *mat64.Dense
+	}{
+		{
+			data: mat64.NewDense(3, 3, []float64{
+				1, 2, 3,
+				3, 4, 5,
+				5, 6, 7,
+			}),
+			weights: nil,
+			ans: mat64.NewDense(3, 3, []float64{
+				1, 1, 1,
+				1, 1, 1,
+				1, 1, 1,
+			}),
+		},
+	} {
+		// Make a copy of the data to check that it isn't changing.
+		r := test.data.RawMatrix()
+		d := make([]float64, len(r.Data))
+		copy(d, r.Data)
+
+		w := make([]float64, len(test.weights))
+		if test.weights != nil {
+			copy(w, test.weights)
+		}
+		c := CorrelationMatrix(nil, test.data, test.weights)
+		if !c.Equals(test.ans) {
+			t.Errorf("%d: expected corr %v, found %v", i, test.ans, c)
+		}
+		if !floats.Equal(d, r.Data) {
+			t.Errorf("%d: data was modified during execution", i)
+		}
+		if !floats.Equal(w, test.weights) {
+			t.Errorf("%d: weights was modified during execution", i)
+		}
+
+		// compare with call to Covariance
+		_, cols := c.Dims()
+		for ci := 0; ci < cols; ci++ {
+			for cj := 0; cj < cols; cj++ {
+				x := test.data.Col(nil, ci)
+				y := test.data.Col(nil, cj)
+				corr := Correlation(x, y, test.weights)
+				if math.Abs(corr-c.At(ci, cj)) > 1e-14 {
+					t.Errorf("CorrMat does not match at (%v, %v). Want %v, got %v.", ci, cj, corr, c.At(ci, cj))
+				}
+			}
+		}
+
+	}
+	if !Panics(func() { CorrelationMatrix(nil, mat64.NewDense(5, 2, nil), []float64{}) }) {
+		t.Errorf("CorrelationMatrix did not panic with weight size mismatch")
+	}
+	if !Panics(func() { CorrelationMatrix(mat64.NewDense(1, 1, nil), mat64.NewDense(5, 2, nil), nil) }) {
+		t.Errorf("CorrelationMatrix did not panic with preallocation size mismatch")
+	}
+	if !Panics(func() { CorrelationMatrix(nil, mat64.NewDense(2, 2, []float64{1, 2, 3, 4}), []float64{1, -1}) }) {
+		t.Errorf("CorrelationMatrix did not panic with negative weights")
+	}
+}
+
 // benchmarks
 
 func randMat(r, c int) mat64.Matrix {