internal/stats/udist_test.go

// Copyright 2015 The Go 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 stats

import (
	"fmt"
	"math"
	"testing"
)

func aeqTable(a, b [][]float64) bool {
	if len(a) != len(b) {
		return false
	}
	for i := range a {
		if len(a[i]) != len(b[i]) {
			return false
		}
		for j := range a[i] {
			// "%f" precision
			if math.Abs(a[i][j]-b[i][j]) >= 0.000001 {
				return false
			}
		}
	}
	return true
}

// U distribution for N=3 up to U=5.
var udist3 = [][]float64{
	//    m=1         2         3
	{0.250000, 0.100000, 0.050000}, // U=0
	{0.500000, 0.200000, 0.100000}, // U=1
	{0.750000, 0.400000, 0.200000}, // U=2
	{1.000000, 0.600000, 0.350000}, // U=3
	{1.000000, 0.800000, 0.500000}, // U=4
	{1.000000, 0.900000, 0.650000}, // U=5
}

// U distribution for N=5 up to U=5.
var udist5 = [][]float64{
	//    m=1         2         3         4         5
	{0.166667, 0.047619, 0.017857, 0.007937, 0.003968}, // U=0
	{0.333333, 0.095238, 0.035714, 0.015873, 0.007937}, // U=1
	{0.500000, 0.190476, 0.071429, 0.031746, 0.015873}, // U=2
	{0.666667, 0.285714, 0.125000, 0.055556, 0.027778}, // U=3
	{0.833333, 0.428571, 0.196429, 0.095238, 0.047619}, // U=4
	{1.000000, 0.571429, 0.285714, 0.142857, 0.075397}, // U=5
}

func TestUDist(t *testing.T) {
	makeTable := func(n int) [][]float64 {
		out := make([][]float64, 6)
		for U := 0; U < 6; U++ {
			out[U] = make([]float64, n)
			for m := 1; m <= n; m++ {
				out[U][m-1] = UDist{N1: m, N2: n}.CDF(float64(U))
			}
		}
		return out
	}
	fmtTable := func(a [][]float64) string {
		out := fmt.Sprintf("%8s", "m=")
		for m := 1; m <= len(a[0]); m++ {
			out += fmt.Sprintf("%9d", m)
		}
		out += "\n"

		for U, row := range a {
			out += fmt.Sprintf("U=%-6d", U)
			for m := 1; m <= len(a[0]); m++ {
				out += fmt.Sprintf(" %f", row[m-1])
			}
			out += "\n"
		}
		return out
	}

	// Compare against tables given in Mann, Whitney (1947).
	got3 := makeTable(3)
	if !aeqTable(got3, udist3) {
		t.Errorf("For n=3, want:\n%sgot:\n%s", fmtTable(udist3), fmtTable(got3))
	}

	got5 := makeTable(5)
	if !aeqTable(got5, udist5) {
		t.Errorf("For n=5, want:\n%sgot:\n%s", fmtTable(udist5), fmtTable(got5))
	}
}

func BenchmarkUDist(b *testing.B) {
	for i := 0; i < b.N; i++ {
		// R uses the exact distribution up to N=50.
		// N*M/2=1250 is the hardest point to get the CDF for.
		UDist{N1: 50, N2: 50}.CDF(1250)
	}
}

func TestUDistTies(t *testing.T) {
	makeTable := func(m, N int, t []int, minx, maxx float64) [][]float64 {
		out := [][]float64{}
		dist := UDist{N1: m, N2: N - m, T: t}
		for x := minx; x <= maxx; x += 0.5 {
			// Convert x from uQt' to uQv'.
			U := x - float64(m*m)/2
			P := dist.CDF(U)
			if len(out) == 0 || !aeq(out[len(out)-1][1], P) {
				out = append(out, []float64{x, P})
			}
		}
		return out
	}
	fmtTable := func(table [][]float64) string {
		out := ""
		for _, row := range table {
			out += fmt.Sprintf("%5.1f %f\n", row[0], row[1])
		}
		return out
	}

	// Compare against Table 1 from Klotz (1966).
	got := makeTable(5, 10, []int{1, 1, 2, 1, 1, 2, 1, 1}, 12.5, 19.5)
	want := [][]float64{
		{12.5, 0.003968}, {13.5, 0.007937},
		{15.0, 0.023810}, {16.5, 0.047619},
		{17.5, 0.071429}, {18.0, 0.087302},
		{19.0, 0.134921}, {19.5, 0.138889},
	}
	if !aeqTable(got, want) {
		t.Errorf("Want:\n%sgot:\n%s", fmtTable(want), fmtTable(got))
	}

	got = makeTable(10, 21, []int{6, 5, 4, 3, 2, 1}, 52, 87)
	want = [][]float64{
		{52.0, 0.000014}, {56.5, 0.000128},
		{57.5, 0.000145}, {60.0, 0.000230},
		{61.0, 0.000400}, {62.0, 0.000740},
		{62.5, 0.000797}, {64.0, 0.000825},
		{64.5, 0.001165}, {65.5, 0.001477},
		{66.5, 0.002498}, {67.0, 0.002725},
		{67.5, 0.002895}, {68.0, 0.003150},
		{68.5, 0.003263}, {69.0, 0.003518},
		{69.5, 0.003603}, {70.0, 0.005648},
		{70.5, 0.005818}, {71.0, 0.006626},
		{71.5, 0.006796}, {72.0, 0.008157},
		{72.5, 0.009688}, {73.0, 0.009801},
		{73.5, 0.010430}, {74.0, 0.011111},
		{74.5, 0.014230}, {75.0, 0.014612},
		{75.5, 0.017249}, {76.0, 0.018307},
		{76.5, 0.020178}, {77.0, 0.022270},
		{77.5, 0.023189}, {78.0, 0.026931},
		{78.5, 0.028207}, {79.0, 0.029979},
		{79.5, 0.030931}, {80.0, 0.038969},
		{80.5, 0.043063}, {81.0, 0.044262},
		{81.5, 0.046389}, {82.0, 0.049581},
		{82.5, 0.056300}, {83.0, 0.058027},
		{83.5, 0.063669}, {84.0, 0.067454},
		{84.5, 0.074122}, {85.0, 0.077425},
		{85.5, 0.083498}, {86.0, 0.094079},
		{86.5, 0.096693}, {87.0, 0.101132},
	}
	if !aeqTable(got, want) {
		t.Errorf("Want:\n%sgot:\n%s", fmtTable(want), fmtTable(got))
	}

	got = makeTable(8, 16, []int{2, 2, 2, 2, 2, 2, 2, 2}, 32, 54)
	want = [][]float64{
		{32.0, 0.000078}, {34.0, 0.000389},
		{36.0, 0.001088}, {38.0, 0.002642},
		{40.0, 0.005905}, {42.0, 0.011500},
		{44.0, 0.021057}, {46.0, 0.035664},
		{48.0, 0.057187}, {50.0, 0.086713},
		{52.0, 0.126263}, {54.0, 0.175369},
	}
	if !aeqTable(got, want) {
		t.Errorf("Want:\n%sgot:\n%s", fmtTable(want), fmtTable(got))
	}

	// Check remaining tables from Klotz against the reference
	// implementation.
	checkRef := func(n1 int, tie []int) {
		wantPMF1, wantCDF1 := udistRef(n1, tie)

		dist := UDist{N1: n1, N2: sumint(tie) - n1, T: tie}
		gotPMF, wantPMF := [][]float64{}, [][]float64{}
		gotCDF, wantCDF := [][]float64{}, [][]float64{}
		N := sumint(tie)
		for U := 0.0; U <= float64(n1*(N-n1)); U += 0.5 {
			gotPMF = append(gotPMF, []float64{U, dist.PMF(U)})
			gotCDF = append(gotCDF, []float64{U, dist.CDF(U)})
			wantPMF = append(wantPMF, []float64{U, wantPMF1[int(U*2)]})
			wantCDF = append(wantCDF, []float64{U, wantCDF1[int(U*2)]})
		}
		if !aeqTable(wantPMF, gotPMF) {
			t.Errorf("For PMF of n1=%v, t=%v, want:\n%sgot:\n%s", n1, tie, fmtTable(wantPMF), fmtTable(gotPMF))
		}
		if !aeqTable(wantCDF, gotCDF) {
			t.Errorf("For CDF of n1=%v, t=%v, want:\n%sgot:\n%s", n1, tie, fmtTable(wantCDF), fmtTable(gotCDF))
		}
	}
	checkRef(5, []int{1, 1, 2, 1, 1, 2, 1, 1})
	checkRef(5, []int{1, 1, 2, 1, 1, 1, 2, 1})
	checkRef(5, []int{1, 3, 1, 2, 1, 1, 1})
	checkRef(8, []int{1, 2, 1, 1, 1, 1, 2, 2, 1, 2})
	checkRef(12, []int{3, 3, 4, 3, 4, 5})
	checkRef(10, []int{1, 2, 3, 4, 5, 6})
}

func BenchmarkUDistTies(b *testing.B) {
	// Worst case: just one tie.
	n := 20
	t := make([]int, 2*n-1)
	for i := range t {
		t[i] = 1
	}
	t[0] = 2

	for i := 0; i < b.N; i++ {
		UDist{N1: n, N2: n, T: t}.CDF(float64(n*n) / 2)
	}
}

func XTestPrintUmemo(t *testing.T) {
	// Reproduce table from Cheung, Klotz.
	ties := []int{4, 5, 3, 4, 6}
	printUmemo(makeUmemo(80, 10, ties), ties)
}

// udistRef computes the PMF and CDF of the U distribution for two
// samples of sizes n1 and sum(t)-n1 with tie vector t. The returned
// pmf and cdf are indexed by 2*U.
//
// This uses the "graphical method" of Klotz (1966). It is very slow
// (Θ(∏ (t[i]+1)) = Ω(2^|t|)), but very correct, and hence useful as a
// reference for testing faster implementations.
func udistRef(n1 int, t []int) (pmf, cdf []float64) {
	// Enumerate all u vectors for which 0 <= u_i <= t_i. Count
	// the number of permutations of two samples of sizes n1 and
	// sum(t)-n1 with tie vector t and accumulate these counts by
	// their U statistics in count[2*U].
	counts := make([]int, 1+2*n1*(sumint(t)-n1))

	u := make([]int, len(t))
	u[0] = -1 // Get enumeration started.
enumu:
	for {
		// Compute the next u vector.
		u[0]++
		for i := 0; i < len(u) && u[i] > t[i]; i++ {
			if i == len(u)-1 {
				// All u vectors have been enumerated.
				break enumu
			}
			// Carry.
			u[i+1]++
			u[i] = 0
		}

		// Is this a legal u vector?
		if sumint(u) != n1 {
			// Klotz (1966) has a method for directly
			// enumerating legal u vectors, but the point
			// of this is to be correct, not fast.
			continue
		}

		// Compute 2*U statistic for this u vector.
		twoU, vsum := 0, 0
		for i, u_i := range u {
			v_i := t[i] - u_i
			// U = U + vsum*u_i + u_i*v_i/2
			twoU += 2*vsum*u_i + u_i*v_i
			vsum += v_i
		}

		// Compute Π choose(t_i, u_i). This is the number of
		// ways of permuting the input sample under u.
		prod := 1
		for i, u_i := range u {
			prod *= int(mathChoose(t[i], u_i) + 0.5)
		}

		// Accumulate the permutations on this u path.
		counts[twoU] += prod

		if false {
			// Print a table in the form of Klotz's
			// "direct enumeration" example.
			//
			// Convert 2U = 2UQV' to UQt' used in Klotz
			// examples.
			UQt := float64(twoU)/2 + float64(n1*n1)/2
			fmt.Printf("%+v %f %-2d\n", u, UQt, prod)
		}
	}

	// Convert counts into probabilities for PMF and CDF.
	pmf = make([]float64, len(counts))
	cdf = make([]float64, len(counts))
	total := int(mathChoose(sumint(t), n1) + 0.5)
	for i, count := range counts {
		pmf[i] = float64(count) / float64(total)
		if i > 0 {
			cdf[i] = cdf[i-1]
		}
		cdf[i] += pmf[i]
	}
	return
}

// printUmemo prints the output of makeUmemo for debugging.
func printUmemo(A []map[ukey]float64, t []int) {
	fmt.Printf("K\tn1\t2*U\tpr\n")
	for K := len(A) - 1; K >= 0; K-- {
		for i, pr := range A[K] {
			_, ref := udistRef(i.n1, t[:K])
			fmt.Printf("%v\t%v\t%v\t%v\t%v\n", K, i.n1, i.twoU, pr, ref[i.twoU])
		}
	}
}