Skip to content

Commit

Permalink
Merge d9e9699 into a8fc114
Browse files Browse the repository at this point in the history
  • Loading branch information
@kortschak
kortschak committed Jun 17, 2017
2 parents a8fc114 + d9e9699 commit 1d4126e
Show file tree
Hide file tree
Showing 4 changed files with 459 additions and 0 deletions.
131 changes: 131 additions & 0 deletions stat/spatial/spatial.go
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
// Copyright ©2017 The gonum 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 spatial // import "gonum.org/v1/gonum/stat/spatial"

import (
"math"

"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/stat"
)

// 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.
//
// G^*_i = num_i / den_i
//
// num_i = \sum_j (w_{ij} x_j) - \bar X \sum_j w_{ij}
// den_i = S \sqrt((n \sum_j w_{ij}^2 - (\sum_j w_{ij})^2)/(n - 1))
// \bar X = (\sum_j x_j) / n
// S = \sqrt((\sum_j x_j^2)/n - (\bar X)^2)
//
// GetisOrdGStar will panic if locality is not a square matrix with dimensions the
// same as the length of data or if i is not a valid index into data.
//
// See doi.org/10.1111%2Fj.1538-4632.1995.tb00912.x.
//
// Weighted Getis-Ord G*i is not currently implemented and GetisOrdGStar will
// panic if weights is not nil.
func GetisOrdGStar(i int, data, weights []float64, locality mat.Matrix) float64 {
if weights != nil {
panic("spatial: weighted data not yet implemented")
}
r, c := locality.Dims()
if r != len(data) || c != len(data) {
panic("spatial: data length mismatch")
}

n := float64(len(data))
mean := stat.Mean(data, weights)
var ss, dwd, dww, sw float64
for j, v := range data {
ss += v * v
w := locality.At(i, j)
sw += w
dwd += w * v
dww += w * w
}
s := math.Sqrt(ss/n - mean*mean)

return (dwd - mean*sw) / (s * math.Sqrt((n*dww-sw*sw)/(n-1)))
}

// GlobalMoransI calculates Global Moran's I calculation of spatial autocorrelation
// for the given data using the provided locality matrix. GlobalMoransI returns
// Moran's I, Var(I) and the z-score associated with those values.
// GlobalMoransI will panic if locality is not a square matrix with dimensions the
// same as the length of data.
//
// See https://doi.org/10.1111%2Fj.1538-4632.2007.00708.x.
//
// Weighted Global Moran's I is not currently implemented and GlobalMoransI will
// panic if weights is not nil.
func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float64) {
if weights != nil {
panic("spatial: weighted data not yet implemented")
}
if r, c := locality.Dims(); r != len(data) || c != len(data) {
panic("spatial: data length mismatch")
}
mean := stat.Mean(data, nil)

// Calculate Moran's I for the data.
var num, den float64
for i, xi := range data {
zi := xi - mean
den += zi * zi
for j, xj := range data {
zj := xj - mean
num += locality.At(i, j) * zi * zj
}
}
i = (float64(len(data)) * num) / (mat.Sum(locality) * den)

// Calculate Moran's E(I) for the data.
e := -1 / float64(len(data)-1)

// Calculate Moran's Var(I) for the data.
// http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm
// http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-global-morans-i-additional-math.htm
var s0, s1, s2 float64
var var2, var4 float64
for i, v := range data {
v -= mean
v *= v
var2 += v
var4 += v * v

var p2 float64
for j := range data {
wij := locality.At(i, j)
wji := locality.At(j, i)

s0 += wij

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

p2 += wij + wji
}
s2 += p2 * p2
}
s1 *= 0.5

n := float64(len(data))
a := n * ((n*n-3*n+3)*s1 - n*s2 + 3*s0*s0)
c := (n - 1) * (n - 2) * (n - 3) * s0 * s0
d := var4 / (var2 * var2)
b := d * ((n*n-n)*s1 - 2*n*s2 + 6*s0*s0)

v = (a-b)/c - e*e

// Calculate z-score associated with Moran's I for the data.
z = (i - e) / math.Sqrt(v)

return i, v, z
}
71 changes: 71 additions & 0 deletions stat/spatial/spatial_areal_example_test.go
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
// Copyright ©2017 The gonum 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 spatial_test

import (
"fmt"
"math"

"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/stat/spatial"
)

type Euclid struct{ x, y int }

func (e Euclid) Dims() (r, c int) { return e.x * e.y, e.x * e.y }
func (e Euclid) At(i, j int) float64 {
d := e.x * e.y
if i < 0 || d <= i || j < 0 || d <= j {
panic("bounds error")
}
if i == j {
return 0
}
x := float64(j%e.x - i%e.x)
y := float64(j/e.x - i/e.x)
return 1 / math.Hypot(x, y)
}
func (e Euclid) T() mat.Matrix { return mat.Transpose{e} }

func ExampleGlobalMoransI_2() {
locality := Euclid{10, 10}

data1 := []float64{
1, 0, 0, 1, 0, 0, 1, 0, 0, 0,
0, 1, 1, 0, 0, 1, 0, 0, 0, 0,
1, 0, 0, 1, 0, 0, 0, 0, 1, 0,
0, 0, 1, 0, 1, 0, 1, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 1, 0, 1, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 1, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 1, 0,
}
i1, _, z1 := spatial.GlobalMoransI(data1, nil, locality)

data2 := []float64{
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
}
i2, _, z2 := spatial.GlobalMoransI(data2, nil, locality)

fmt.Printf("%v sparse points Moran's I=%.4v z-score=%.4v\n", floats.Sum(data1), i1, z1)
fmt.Printf("%v clustered points Moran's I=%.4v z-score=%.4v\n", floats.Sum(data2), i2, z2)

// Output:
//
// 24 sparse points Moran's I=-0.02999 z-score=-1.913
// 24 clustered points Moran's I=0.09922 z-score=10.52
}
69 changes: 69 additions & 0 deletions stat/spatial/spatial_example_test.go
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
// Copyright ©2017 The gonum 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 spatial_test

import (
"fmt"

"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/stat/spatial"
)

func ExampleGlobalMoransI_1() {
data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}
locality := mat.NewDense(10, 10, []float64{
0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 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}
locality := mat.NewDense(10, 10, []float64{
1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
})

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
}

0 comments on commit 1d4126e

Please sign in to comment.