Merge 01acf26 into e15809c

stat/spatial.go

@@ -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())
+}

stat/spatial_areal_example_test.go

@@ -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
+}

stat/spatial_example_test.go

@@ -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
+}