stat/spatial: use mat.RowNonZeroDoer
MODIFY
stat/spatial/spatial.go
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,12 +42,20 @@ 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
+ 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)
return (dwd - mean*sw) / (s * math.Sqrt((n*dww-sw*sw)/(n-1)))
@@ -73,11 +80,20 @@ 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
+ 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
@@ -85,6 +101,7 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
num += w * zi * zj
}
}
+ }
i = (float64(len(data)) / sum) * (num / den)
// Calculate Moran's E(I) for the data.
@@ -102,6 +119,18 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
var4 += v * v
var p2 float64
+ if isDoer {
+ doer.DoRowNonZero(i, func(i, j int, wij float64) {
+ wji := locality.At(j, i)
+
+ s0 += wij
+
+ v := wij + wji
+ s1 += v * v
+
+ p2 += v
+ })
+ } else {
for j := range data {
wij := locality.At(i, j)
wji := locality.At(j, i)
@@ -113,6 +142,7 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
p2 += v
}
+ }
s2 += p2 * p2
}
s1 *= 0.5

MODIFY
stat/spatial/spatial_example_test.go
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
+}

MODIFY
stat/spatial/spatial_test.go
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) {

stat/spatial: use mat.RowNonZeroDoer