1// Copyright ©2016 The Gonum Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package stat_test
6
7import (
8	"testing"
9
10	"gonum.org/v1/gonum/floats"
11	"gonum.org/v1/gonum/mat"
12	"gonum.org/v1/gonum/stat"
13)
14
15func TestCanonicalCorrelations(t *testing.T) {
16tests:
17	for i, test := range []struct {
18		xdata     mat.Matrix
19		ydata     mat.Matrix
20		weights   []float64
21		wantCorrs []float64
22		wantpVecs *mat.Dense
23		wantqVecs *mat.Dense
24		wantphiVs *mat.Dense
25		wantpsiVs *mat.Dense
26		epsilon   float64
27	}{
28		// Test results verified using R.
29		{ // Truncated iris data, Sepal vs Petal measurements.
30			xdata: mat.NewDense(10, 2, []float64{
31				5.1, 3.5,
32				4.9, 3.0,
33				4.7, 3.2,
34				4.6, 3.1,
35				5.0, 3.6,
36				5.4, 3.9,
37				4.6, 3.4,
38				5.0, 3.4,
39				4.4, 2.9,
40				4.9, 3.1,
41			}),
42			ydata: mat.NewDense(10, 2, []float64{
43				1.4, 0.2,
44				1.4, 0.2,
45				1.3, 0.2,
46				1.5, 0.2,
47				1.4, 0.2,
48				1.7, 0.4,
49				1.4, 0.3,
50				1.5, 0.2,
51				1.4, 0.2,
52				1.5, 0.1,
53			}),
54			wantCorrs: []float64{0.7250624174504773, 0.5547679185730191},
55			wantpVecs: mat.NewDense(2, 2, []float64{
56				0.0765914610875867, 0.9970625597666721,
57				0.9970625597666721, -0.0765914610875868,
58			}),
59			wantqVecs: mat.NewDense(2, 2, []float64{
60				0.3075184850910837, 0.9515421069649439,
61				0.9515421069649439, -0.3075184850910837,
62			}),
63			wantphiVs: mat.NewDense(2, 2, []float64{
64				-1.9794877596804641, 5.2016325219025124,
65				4.5211829944066553, -2.7263663170835697,
66			}),
67			wantpsiVs: mat.NewDense(2, 2, []float64{
68				-0.0613084818030103, 10.8514169865438941,
69				12.7209032660734298, -7.6793888180353775,
70			}),
71			epsilon: 1e-12,
72		},
73		// Test results compared to those results presented in examples by
74		// Koch, Inge. Analysis of multivariate and high-dimensional data.
75		// Vol. 32. Cambridge University Press, 2013. ISBN: 9780521887939
76		{ // ASA Car Exposition Data of Ramos and Donoho (1983)
77			// Displacement, Horsepower, Weight
78			xdata: carData.Slice(0, 392, 0, 3),
79			// Acceleration, MPG
80			ydata:     carData.Slice(0, 392, 3, 5),
81			wantCorrs: []float64{0.8782187384352336, 0.6328187219216761},
82			wantpVecs: mat.NewDense(3, 2, []float64{
83				0.3218296374829181, 0.3947540257657075,
84				0.4162807660635797, 0.7573719053303306,
85				0.8503740401982725, -0.5201509936144236,
86			}),
87			wantqVecs: mat.NewDense(2, 2, []float64{
88				-0.5161984172278830, -0.8564690269072364,
89				-0.8564690269072364, 0.5161984172278830,
90			}),
91			wantphiVs: mat.NewDense(3, 2, []float64{
92				0.0025033152994308, 0.0047795464118615,
93				0.0201923608080173, 0.0409150208725958,
94				-0.0000247374128745, -0.0026766435161875,
95			}),
96			wantpsiVs: mat.NewDense(2, 2, []float64{
97				-0.1666196759760772, -0.3637393866139658,
98				-0.0915512109649727, 0.1077863777929168,
99			}),
100			epsilon: 1e-12,
101		},
102		// Test results compared to those results presented in examples by
103		// Koch, Inge. Analysis of multivariate and high-dimensional data.
104		// Vol. 32. Cambridge University Press, 2013. ISBN: 9780521887939
105		{ // Boston Housing Data of Harrison and Rubinfeld (1978)
106			// Per capita crime rate by town,
107			// Proportion of non-retail business acres per town,
108			// Nitric oxide concentration (parts per 10 million),
109			// Weighted distances to Boston employment centres,
110			// Index of accessibility to radial highways,
111			// Pupil-teacher ratio by town, Proportion of blacks by town
112			xdata: bostonData.Slice(0, 506, 0, 7),
113			// Average number of rooms per dwelling,
114			// Proportion of owner-occupied units built prior to 1940,
115			// Full-value property-tax rate per $10000,
116			// Median value of owner-occupied homes in $1000s
117			ydata:     bostonData.Slice(0, 506, 7, 11),
118			wantCorrs: []float64{0.9451239443886021, 0.6786622733370654, 0.5714338361583764, 0.2009739704710440},
119			wantpVecs: mat.NewDense(7, 4, []float64{
120				-0.2574391924541903, 0.0158477516621194, 0.2122169934631024, -0.0945733803894706,
121				-0.4836594430018478, 0.3837101908138468, 0.1474448317415911, 0.6597324886718275,
122				-0.0800776365873296, 0.3493556742809252, 0.3287336458109373, -0.2862040444334655,
123				0.1277586360386374, -0.7337427663667596, 0.4851134819037011, 0.2247964865970192,
124				-0.6969432006136684, -0.4341748776002893, -0.3602872887636357, 0.0290661608626292,
125				-0.0990903250057199, 0.0503411215453873, 0.6384330631742202, 0.1022367136218303,
126				0.4260459963765036, 0.0323334351308141, -0.2289527516030810, 0.6419232947608805,
127			}),
128			wantqVecs: mat.NewDense(4, 4, []float64{
129				0.0181660502363264, -0.1583489460479038, -0.0066723577642883, -0.9871935400650649,
130				-0.2347699045986119, 0.9483314614936594, -0.1462420505631345, -0.1554470767919033,
131				-0.9700704038477141, -0.2406071741000039, -0.0251838984227037, 0.0209134074358349,
132				0.0593000682318482, -0.1330460003097728, -0.9889057151969489, 0.0291161494720761,
133			}),
134			wantphiVs: mat.NewDense(7, 4, []float64{
135				-0.0027462234108197, 0.0093444513500898, 0.0489643932714296, -0.0154967189805819,
136				-0.0428564455279537, -0.0241708702119420, 0.0360723472093996, 0.1838983230588095,
137				-1.2248435648802380, 5.6030921364723980, 5.8094144583797025, -4.7926812190419676,
138				-0.0043684825094649, -0.3424101164977618, 0.4469961215717917, 0.1150161814353696,
139				-0.0741534069521954, -0.1193135794923700, -0.1115518305471460, 0.0021638758323088,
140				-0.0233270323101624, 0.1046330818178399, 0.3853045975077387, -0.0160927870102877,
141				0.0001293051387859, 0.0004540746921446, -0.0030296315865440, 0.0081895477974654,
142			}),
143			wantpsiVs: mat.NewDense(4, 4, []float64{
144				0.0301593362017375, -0.3002219289647127, 0.0878217377593682, -1.9583226531517062,
145				-0.0065483104073892, 0.0392212086716247, -0.0117570776209991, -0.0061113064481860,
146				-0.0052075523350125, -0.0045770200452960, -0.0022762313289592, 0.0008441873006821,
147				0.0020111735096327, 0.0037352799829930, -0.1292578071621794, 0.1037709056329765,
148			}),
149			epsilon: 1e-12,
150		},
151	} {
152		var cc stat.CC
153		var corrs []float64
154		var pVecs, qVecs mat.Dense
155		var phiVs, psiVs mat.Dense
156		for j := 0; j < 2; j++ {
157			err := cc.CanonicalCorrelations(test.xdata, test.ydata, test.weights)
158			if err != nil {
159				t.Errorf("%d use %d: unexpected error: %v", i, j, err)
160				continue tests
161			}
162
163			corrs = cc.CorrsTo(corrs)
164			cc.LeftTo(&pVecs, true)
165			cc.RightTo(&qVecs, true)
166			cc.LeftTo(&phiVs, false)
167			cc.RightTo(&psiVs, false)
168
169			if !floats.EqualApprox(corrs, test.wantCorrs, test.epsilon) {
170				t.Errorf("%d use %d: unexpected variance result got:%v, want:%v",
171					i, j, corrs, test.wantCorrs)
172			}
173			if !mat.EqualApprox(&pVecs, test.wantpVecs, test.epsilon) {
174				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
175					i, j, mat.Formatted(&pVecs), mat.Formatted(test.wantpVecs))
176			}
177			if !mat.EqualApprox(&qVecs, test.wantqVecs, test.epsilon) {
178				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
179					i, j, mat.Formatted(&qVecs), mat.Formatted(test.wantqVecs))
180			}
181			if !mat.EqualApprox(&phiVs, test.wantphiVs, test.epsilon) {
182				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
183					i, j, mat.Formatted(&phiVs), mat.Formatted(test.wantphiVs))
184			}
185			if !mat.EqualApprox(&psiVs, test.wantpsiVs, test.epsilon) {
186				t.Errorf("%d use %d: unexpected CCA result got:\n%v\nwant:\n%v",
187					i, j, mat.Formatted(&psiVs), mat.Formatted(test.wantpsiVs))
188			}
189		}
190	}
191}
192