-
Notifications
You must be signed in to change notification settings - Fork 535
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
stat: add spatial statistic measures
- Loading branch information
Showing
4 changed files
with
481 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,148 @@ | ||
// 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 stat | ||
|
||
import ( | ||
"math" | ||
|
||
"gonum.org/v1/gonum/matrix/mat64" | ||
) | ||
|
||
// GetisOrd performs Local Getis-Ord G*i statistic calculation. | ||
type GetisOrd struct { | ||
data *mat64.Vector | ||
locality *mat64.Dense | ||
|
||
mean, s float64 | ||
} | ||
|
||
// NewGetisOrd returns a GetisOrd based on the provided data and locality matrix. | ||
// NewGetisOrd will panic if locality is not a square matrix with dimensions the | ||
// same as the length of data. | ||
func NewGetisOrd(data []float64, locality mat64.Matrix) GetisOrd { | ||
r, c := locality.Dims() | ||
if r != len(data) || c != len(data) { | ||
panic("stat: data length mismatch") | ||
} | ||
|
||
var g GetisOrd | ||
d := make([]float64, len(data)) | ||
copy(d, data) | ||
g.data = mat64.NewVector(len(data), d) | ||
|
||
g.locality = mat64.DenseCopyOf(locality) | ||
g.mean = Mean(data, nil) | ||
var ss float64 | ||
for _, v := range data { | ||
ss += v * v | ||
} | ||
g.s = ss/float64(len(data)) - g.mean*g.mean | ||
|
||
return g | ||
} | ||
|
||
// Len returns the number of data points held by the receiver. | ||
func (g GetisOrd) Len() int { return g.data.Len() } | ||
|
||
// Gstar returns the Local Getis-Ord G* statistic for element i. The | ||
// returned value is a z-score. | ||
func (g GetisOrd) Gstar(i int) float64 { | ||
wi := g.locality.RowView(i) | ||
ws := g.mean * mat64.Sum(wi) | ||
n := float64(g.data.Len()) | ||
num := mat64.Dot(wi, g.data) - ws | ||
den := g.s * math.Sqrt((n*mat64.Dot(wi, wi)-ws*ws)/(n-1)) | ||
return num / den | ||
} | ||
|
||
// Moran performs Global Moran's I calculation of spatial autocorrelation. | ||
// | ||
// See https://en.wikipedia.org/wiki/Moran%27s_I. | ||
type Moran struct { | ||
data []float64 | ||
mean float64 | ||
locality *mat64.Dense | ||
} | ||
|
||
// NewMoran returns a new Moran based on the provided data and locality matrix. | ||
// NewMoran will panic if locality is not a square matrix with dimensions the | ||
// same as the length of data. | ||
func NewMoran(data []float64, locality mat64.Matrix) Moran { | ||
r, c := locality.Dims() | ||
if r != len(data) || c != len(data) { | ||
panic("stat: data length mismatch") | ||
} | ||
var m Moran | ||
m.data = make([]float64, len(data)) | ||
copy(m.data, data) | ||
m.mean = Mean(data, nil) | ||
m.locality = mat64.DenseCopyOf(locality) | ||
|
||
return m | ||
} | ||
|
||
// I returns Moran's I for the data represented by the receiver. | ||
func (m Moran) I() float64 { | ||
var num, den float64 | ||
for i, xi := range m.data { | ||
zi := xi - m.mean | ||
den += zi * zi | ||
for j, xj := range m.data { | ||
zj := xj - m.mean | ||
num += m.locality.At(i, j) * zi * zj | ||
} | ||
} | ||
return (float64(len(m.data)) * num) / (mat64.Sum(m.locality) * den) | ||
} | ||
|
||
// E returns Moran's E(I) for the data represented by the receiver. | ||
func (m Moran) E() float64 { | ||
return -1 / float64(len(m.data)-1) | ||
} | ||
|
||
// Var returns Moran's Var(I) for the data represented by the receiver. | ||
func (m Moran) Var() float64 { | ||
// From http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm | ||
// and 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 m.data { | ||
v -= m.mean | ||
v *= v | ||
var2 += v | ||
var4 += v * v | ||
|
||
var p2 float64 | ||
for j := range m.data { | ||
wij := m.locality.At(i, j) | ||
wji := m.locality.At(j, i) | ||
|
||
s0 += wij | ||
|
||
v := wij + wji | ||
s1 += v * v | ||
|
||
p2 += wij + wji | ||
} | ||
s2 += p2 * p2 | ||
} | ||
s1 *= 0.5 | ||
|
||
n := float64(len(m.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) | ||
|
||
e := m.E() | ||
|
||
return (a-b)/c - e*e | ||
} | ||
|
||
// Z returns the z-score associated with Moran's I for the data represented by the receiver. | ||
func (m Moran) Z() float64 { | ||
return (m.I() - m.E()) / math.Sqrt(m.Var()) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 stat_test | ||
|
||
import ( | ||
"fmt" | ||
"math" | ||
|
||
"gonum.org/v1/gonum/floats" | ||
"gonum.org/v1/gonum/matrix/mat64" | ||
"gonum.org/v1/gonum/stat" | ||
) | ||
|
||
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() mat64.Matrix { return mat64.Transpose{e} } | ||
|
||
func ExampleMoran_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, | ||
} | ||
m1 := stat.NewMoran(data1, 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, | ||
} | ||
m2 := stat.NewMoran(data2, locality) | ||
|
||
fmt.Printf("%v sparse points Moran's I=%.4v z-score=%.4v\n", floats.Sum(data1), m1.I(), m1.Z()) | ||
fmt.Printf("%v clustered points Moran's I=%.4v z-score=%.4v\n", floats.Sum(data2), m2.I(), m2.Z()) | ||
|
||
// 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 | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 stat_test | ||
|
||
import ( | ||
"fmt" | ||
|
||
"gonum.org/v1/gonum/matrix/mat64" | ||
"gonum.org/v1/gonum/stat" | ||
) | ||
|
||
func ExampleMoran_1() { | ||
data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0} | ||
locality := mat64.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, | ||
}) | ||
|
||
m := stat.NewMoran(data, locality) | ||
|
||
fmt.Printf("Moran's I=%.4v z-score=%.4v\n", m.I(), m.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 := mat64.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, | ||
}) | ||
|
||
g := stat.NewGetisOrd(data, locality) | ||
|
||
for i, v := range data { | ||
fmt.Printf("v=%v G*i=% .4v\n", v, g.Gstar(i)) | ||
} | ||
|
||
// Output: | ||
// | ||
// v=0 G*i=-2.273 | ||
// v=0 G*i=-2.807 | ||
// v=0 G*i=-0.4678 | ||
// v=1 G*i= 1.871 | ||
// v=1 G*i= 4.21 | ||
// v=1 G*i= 1.871 | ||
// v=0 G*i= 1.871 | ||
// v=1 G*i=-0.4678 | ||
// v=0 G*i=-0.4678 | ||
// v=0 G*i=-2.273 | ||
} |
Oops, something went wrong.