diff --git a/stat/spatial/spatial.go b/stat/spatial/spatial.go index eca484dc..a45fe422 100644 --- a/stat/spatial/spatial.go +++ b/stat/spatial/spatial.go @@ -12,7 +12,6 @@ import ( ) // TODO(kortschak): Implement weighted routines. -// TODO(kortschak): Make use of banded matrices when they exist in mat. // GetisOrdGStar returns the Local Getis-Ord G*i statistic for element of of the // weighted data using the provided locality matrix. The returned value is a z-score. @@ -43,11 +42,19 @@ func GetisOrdGStar(i int, data, weights []float64, locality mat.Matrix) float64 n := float64(len(data)) mean, std := stat.MeanStdDev(data, weights) var dwd, dww, sw float64 - for j, v := range data { - w := locality.At(i, j) - sw += w - dwd += w * v - dww += w * w + if doer, ok := locality.(mat.RowNonZeroDoer); ok { + doer.DoRowNonZero(i, func(_, j int, w float64) { + sw += w + dwd += w * data[j] + dww += w * w + }) + } else { + for j, v := range data { + w := locality.At(i, j) + sw += w + dwd += w * v + dww += w * w + } } s := std * math.Sqrt((n-1)/n) @@ -73,16 +80,26 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6 } mean := stat.Mean(data, nil) + doer, isDoer := locality.(mat.RowNonZeroDoer) + // Calculate Moran's I for the data. var num, den, sum float64 for i, xi := range data { zi := xi - mean den += zi * zi - for j, xj := range data { - w := locality.At(i, j) - sum += w - zj := xj - mean - num += w * zi * zj + if isDoer { + doer.DoRowNonZero(i, func(_, j int, w float64) { + sum += w + zj := data[j] - mean + num += w * zi * zj + }) + } else { + for j, xj := range data { + w := locality.At(i, j) + sum += w + zj := xj - mean + num += w * zi * zj + } } } i = (float64(len(data)) / sum) * (num / den) @@ -102,16 +119,29 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6 var4 += v * v var p2 float64 - for j := range data { - wij := locality.At(i, j) - wji := locality.At(j, i) + if isDoer { + doer.DoRowNonZero(i, func(i, j int, wij float64) { + wji := locality.At(j, i) - s0 += wij + s0 += wij - v := wij + wji - s1 += v * v + v := wij + wji + s1 += v * v - p2 += v + p2 += v + }) + } else { + for j := range data { + wij := locality.At(i, j) + wji := locality.At(j, i) + + s0 += wij + + v := wij + wji + s1 += v * v + + p2 += v + } } s2 += p2 * p2 } diff --git a/stat/spatial/spatial_example_test.go b/stat/spatial/spatial_example_test.go index 31a2330a..e5281a7e 100644 --- a/stat/spatial/spatial_example_test.go +++ b/stat/spatial/spatial_example_test.go @@ -38,6 +38,35 @@ func ExampleGlobalMoransI_linear() { // Moran's I=0.1111 z-score=0.6335 } +func ExampleGlobalMoransI_banded() { + data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0} + + // The locality here describes spatial neighbor + // relationships. + // This example uses the band matrix representation + // to improve time and space efficiency. + locality := mat.NewBandDense(10, 10, 1, 1, []float64{ + 0, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 1, + 1, 0, 0, + }) + + i, _, z := spatial.GlobalMoransI(data, nil, locality) + + fmt.Printf("Moran's I=%.4v z-score=%.4v\n", i, z) + + // Output: + // + // Moran's I=0.1111 z-score=0.6335 +} + func ExampleGetisOrdGStar() { data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0} @@ -73,3 +102,41 @@ func ExampleGetisOrdGStar() { // v=0 G*i=-0.2673 // v=0 G*i=-1.225 } + +func ExampleGetisOrd_band() { + data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0} + + // The locality here describes spatial neighbor + // relationships including self. + // This example uses the band matrix representation + // to improve time and space efficiency. + locality := mat.NewBandDense(10, 10, 1, 1, []float64{ + 0, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 1, + 1, 1, 0, + }) + + for i, v := range data { + fmt.Printf("v=%v G*i=% .4v\n", v, spatial.GetisOrdGStar(i, data, nil, locality)) + } + + // Output: + // + // v=0 G*i=-1.225 + // v=0 G*i=-1.604 + // v=0 G*i=-0.2673 + // v=1 G*i= 1.069 + // v=1 G*i= 2.405 + // v=1 G*i= 1.069 + // v=0 G*i= 1.069 + // v=1 G*i=-0.2673 + // v=0 G*i=-0.2673 + // v=0 G*i=-1.225 +} diff --git a/stat/spatial/spatial_test.go b/stat/spatial/spatial_test.go index 1e281974..f91edfdf 100644 --- a/stat/spatial/spatial_test.go +++ b/stat/spatial/spatial_test.go @@ -13,7 +13,7 @@ import ( "gonum.org/v1/gonum/mat" ) -func simpleAdjacency(n, wide int, diag bool) *mat.Dense { +func simpleAdjacency(n, wide int, diag bool) mat.Matrix { m := mat.NewDense(n, n, nil) for i := 0; i < n; i++ { for j := 1; j <= wide; j++ { @@ -30,11 +30,28 @@ func simpleAdjacency(n, wide int, diag bool) *mat.Dense { return m } +func simpleAdjacencyBand(n, wide int, diag bool) mat.Matrix { + m := mat.NewBandDense(n, n, wide, wide, nil) + for i := 0; i < n; i++ { + for j := 1; j <= wide; j++ { + if j > i { + continue + } + m.SetBand(i-j, i, 1) + m.SetBand(i, i-j, 1) + } + if diag { + m.SetBand(i, i, 1) + } + } + return m +} + var spatialTests = []struct { from, to float64 n, wide int fn func(float64, int, *rand.Rand) float64 - locality func(n, wide int, diag bool) *mat.Dense + locality func(n, wide int, diag bool) mat.Matrix // Values for MoranI and z-score are obtained from // an R reference implementation. @@ -46,6 +63,7 @@ var spatialTests = []struct { // of the plotted data. wantSegs int }{ + // Dense matrix locality. { from: 0, to: 1, n: 1000, wide: 1, fn: func(_ float64, _ int, rnd *rand.Rand) float64 { @@ -122,6 +140,84 @@ var spatialTests = []struct { wantZ: -31.559531064275987, wantSegs: 0, }, + + // Band matrix locality. + { + from: 0, to: 1, n: 1000, wide: 1, + fn: func(_ float64, _ int, rnd *rand.Rand) float64 { + return rnd.Float64() + }, + locality: simpleAdjacencyBand, + + wantMoranI: -0.0019631298955953233, + wantZ: -0.03039477405151108, + wantSegs: 0, + }, + { + from: -math.Pi / 2, to: 3 * math.Pi / 2, n: 1000, wide: 1, + fn: func(x float64, _ int, _ *rand.Rand) float64 { + y := math.Sin(x) + if math.Abs(y) > 0.5 { + y *= 1/math.Abs(y) - 1 + } + return y * math.Sin(x*2) + }, + locality: simpleAdjacencyBand, + + wantMoranI: 1.0008149537991464, + wantZ: 31.648547078779092, + wantSegs: 4, + }, + { + from: 0, to: 1, n: 1000, wide: 1, + fn: func(_ float64, _ int, rnd *rand.Rand) float64 { + return rnd.NormFloat64() + }, + locality: simpleAdjacencyBand, + + wantMoranI: 0.031195199553564902, + wantZ: 1.0171161514080056, + wantSegs: 0, + }, + { + from: 0, to: 1, n: 1000, wide: 1, + fn: func(x float64, _ int, rnd *rand.Rand) float64 { + if rnd.Float64() < 0.5 { + return rnd.NormFloat64() + 5 + } + return rnd.NormFloat64() + }, + locality: simpleAdjacencyBand, + + wantMoranI: -0.016245135637562223, + wantZ: -0.48157993864993476, + wantSegs: 0, + }, + { + from: 0, to: 1, n: 1000, wide: 1, + fn: func(x float64, i int, rnd *rand.Rand) float64 { + if i%2 == 0 { + return rnd.NormFloat64() + 5 + } + return rnd.NormFloat64() + }, + locality: simpleAdjacencyBand, + + wantMoranI: -0.8565268969272998, + wantZ: -27.027057520918113, + wantSegs: 0, + }, + { + from: 0, to: 1, n: 1000, wide: 1, + fn: func(_ float64, i int, _ *rand.Rand) float64 { + return float64(i % 2) + }, + locality: simpleAdjacencyBand, + + wantMoranI: -1, + wantZ: -31.559531064275987, + wantSegs: 0, + }, } func TestGetisOrd(t *testing.T) {