Skip to content

Commit

Permalink
stat/spatial: use mat.Banded
Browse files Browse the repository at this point in the history
  • Loading branch information
@kortschak
kortschak committed Jul 7, 2017
1 parent 388406a commit a58dd46
Show file tree
Hide file tree
Showing 3 changed files with 237 additions and 21 deletions.
91 changes: 72 additions & 19 deletions stat/spatial/spatial.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -43,11 +42,21 @@ 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 band, ok := locality.(mat.Banded); ok {
kl, ku := band.Bandwidth()
for j := max(0, i-kl); j < min(c, i+ku+1); j++ {
w := locality.At(i, j)
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)

Expand All @@ -68,21 +77,37 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6
if weights != nil {
panic("spatial: weighted data not yet implemented")
}
if r, c := locality.Dims(); r != len(data) || c != len(data) {
rows, cols := locality.Dims()
if rows != len(data) || cols != len(data) {
panic("spatial: data length mismatch")
}
mean := stat.Mean(data, nil)

var kl, ku int
band, isBand := locality.(mat.Banded)
if isBand {
kl, ku = band.Bandwidth()
}

// 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 isBand {
for j := max(0, i-kl); j < min(cols, i+ku+1); j++ {
w := locality.At(i, j)
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)
Expand All @@ -102,16 +127,30 @@ 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 isBand {
for j := max(0, i-kl); j < min(cols, i+ku+1); j++ {
wij := locality.At(i, j)
wji := locality.At(j, i)

s0 += wij

s0 += wij
v := wij + wji
s1 += v * v

v := wij + wji
s1 += v * v
p2 += v
}
} else {
for j := range data {
wij := locality.At(i, j)
wji := locality.At(j, i)

p2 += v
s0 += wij

v := wij + wji
s1 += v * v

p2 += v
}
}
s2 += p2 * p2
}
Expand All @@ -130,3 +169,17 @@ func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float6

return i, v, z
}

func min(a, b int) int {
if a < b {
return a
}
return b
}

func max(a, b int) int {
if a > b {
return a
}
return b
}
67 changes: 67 additions & 0 deletions stat/spatial/spatial_example_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ExampleGetisOrd() {
data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}

Expand Down Expand Up @@ -73,3 +102,41 @@ func ExampleGetisOrd() {
// 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
}
100 changes: 98 additions & 2 deletions stat/spatial/spatial_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++ {
Expand All @@ -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.
Expand All @@ -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 {
Expand Down Expand Up @@ -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) {
Expand Down

0 comments on commit a58dd46

Please sign in to comment.