mirror of
https://github.com/gonum/gonum.git
synced 2025-10-06 15:47:01 +08:00
stat/spatial: new package for spatial statistic measures
This commit is contained in:
kortschak
132
stat/spatial/spatial.go
Normal file
132
stat/spatial/spatial.go
Normal file
@@ -0,0 +1,132 @@
|
||||
// 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, 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
|
||||
}
|
||||
s := std * math.Sqrt((n-1)/n)
|
||||
|
||||
return (dwd - mean*sw) / (s * math.Sqrt((n*dww-sw*sw)/(n-1)))
|
||||
}
|
||||
|
||||
// GlobalMoransI performs 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, 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
|
||||
}
|
||||
}
|
||||
i = (float64(len(data)) / sum) * (num / 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 += v
|
||||
}
|
||||
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
|
||||
}
|
||||
74
stat/spatial/spatial_areal_example_test.go
Normal file
74
stat/spatial/spatial_areal_example_test.go
Normal file
@@ -0,0 +1,74 @@
|
||||
// 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"
|
||||
)
|
||||
|
||||
// Euclid is a mat.Matrix whose elements refects the Euclidean
|
||||
// distance between a series of unit-separated points strided
|
||||
// to be arranged in an x by y grid.
|
||||
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_areal() {
|
||||
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 scattered 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 scattered points Moran's I=-0.02999 z-score=-1.913
|
||||
// 24 clustered points Moran's I=0.09922 z-score=10.52
|
||||
}
|
||||
75
stat/spatial/spatial_example_test.go
Normal file
75
stat/spatial/spatial_example_test.go
Normal file
@@ -0,0 +1,75 @@
|
||||
// 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_linear() {
|
||||
data := []float64{0, 0, 0, 1, 1, 1, 0, 1, 0, 0}
|
||||
|
||||
// The locality here describes spatial neighbor
|
||||
// relationships.
|
||||
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}
|
||||
|
||||
// The locality here describes spatial neighbor
|
||||
// relationships including self.
|
||||
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
|
||||
}
|
||||
197
stat/spatial/spatial_test.go
Normal file
197
stat/spatial/spatial_test.go
Normal file
@@ -0,0 +1,197 @@
|
||||
// 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 (
|
||||
"math"
|
||||
"math/rand"
|
||||
"testing"
|
||||
|
||||
"gonum.org/v1/gonum/floats"
|
||||
"gonum.org/v1/gonum/mat"
|
||||
)
|
||||
|
||||
func simpleAdjacency(n, wide int, diag bool) *mat.Dense {
|
||||
m := mat.NewDense(n, n, nil)
|
||||
for i := 0; i < n; i++ {
|
||||
for j := 1; j <= wide; j++ {
|
||||
if j > i {
|
||||
continue
|
||||
}
|
||||
m.Set(i-j, i, 1)
|
||||
m.Set(i, i-j, 1)
|
||||
}
|
||||
if diag {
|
||||
m.Set(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
|
||||
|
||||
// Values for MoranI and z-score are obtained from
|
||||
// an R reference implementation.
|
||||
wantMoranI float64
|
||||
wantZ float64
|
||||
|
||||
// The value for expected number of significant
|
||||
// segments is obtained from visual inspection
|
||||
// of the plotted data.
|
||||
wantSegs int
|
||||
}{
|
||||
{
|
||||
from: 0, to: 1, n: 1000, wide: 1,
|
||||
fn: func(_ float64, _ int, rnd *rand.Rand) float64 {
|
||||
return rnd.Float64()
|
||||
},
|
||||
locality: simpleAdjacency,
|
||||
|
||||
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: simpleAdjacency,
|
||||
|
||||
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: simpleAdjacency,
|
||||
|
||||
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: simpleAdjacency,
|
||||
|
||||
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: simpleAdjacency,
|
||||
|
||||
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: simpleAdjacency,
|
||||
|
||||
wantMoranI: -1,
|
||||
wantZ: -31.559531064275987,
|
||||
wantSegs: 0,
|
||||
},
|
||||
}
|
||||
|
||||
func TestGetisOrd(t *testing.T) {
|
||||
for ti, test := range spatialTests {
|
||||
rnd := rand.New(rand.NewSource(1))
|
||||
data := make([]float64, test.n)
|
||||
step := (test.to - test.from) / float64(test.n)
|
||||
for i := range data {
|
||||
data[i] = test.fn(test.from+step*float64(i), i, rnd)
|
||||
}
|
||||
locality := test.locality(test.n, test.wide, true)
|
||||
|
||||
nseg := getisOrdSegments(data, nil, locality)
|
||||
if nseg != test.wantSegs {
|
||||
t.Errorf("unexpected number of significant segments for test %d: got:%d want:%d",
|
||||
ti, nseg, test.wantSegs)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// getisOrdSegments returns the number of contiguously significant G*i segemtns in
|
||||
// data. This allows an intuitive validation of the function in lieu of a reference
|
||||
// implementation.
|
||||
func getisOrdSegments(data, weight []float64, locality mat.Matrix) int {
|
||||
const thresh = 2
|
||||
var nseg int
|
||||
segstart := -1
|
||||
for i := range data {
|
||||
gi := GetisOrdGStar(i, data, weight, locality)
|
||||
if segstart != -1 {
|
||||
if math.Abs(gi) < thresh {
|
||||
// Filter short segments.
|
||||
if i-segstart < 5 {
|
||||
segstart = -1
|
||||
continue
|
||||
}
|
||||
|
||||
segstart = -1
|
||||
nseg++
|
||||
}
|
||||
continue
|
||||
}
|
||||
if math.Abs(gi) >= thresh {
|
||||
segstart = i
|
||||
}
|
||||
}
|
||||
if segstart != -1 && len(data)-segstart >= 5 {
|
||||
nseg++
|
||||
}
|
||||
return nseg
|
||||
}
|
||||
|
||||
func TestGlobalMoransI(t *testing.T) {
|
||||
const tol = 1e-14
|
||||
for ti, test := range spatialTests {
|
||||
rnd := rand.New(rand.NewSource(1))
|
||||
data := make([]float64, test.n)
|
||||
step := (test.to - test.from) / float64(test.n)
|
||||
for i := range data {
|
||||
data[i] = test.fn(test.from+step*float64(i), i, rnd)
|
||||
}
|
||||
locality := test.locality(test.n, test.wide, false)
|
||||
|
||||
gotI, _, gotZ := GlobalMoransI(data, nil, locality)
|
||||
|
||||
if !floats.EqualWithinAbsOrRel(gotI, test.wantMoranI, tol, tol) {
|
||||
t.Errorf("unexpected Moran's I value for test %d: got:%v want:%v", ti, gotI, test.wantMoranI)
|
||||
}
|
||||
if !floats.EqualWithinAbsOrRel(gotZ, test.wantZ, tol, tol) {
|
||||
t.Errorf("unexpected Moran's I z-score for test %d: got:%v want:%v", ti, gotZ, test.wantZ)
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user