diff --git a/stat/spatial.go b/stat/spatial.go
new file mode 100644
index 0000000000..a342845883
--- /dev/null
+++ b/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())
+}
diff --git a/stat/spatial_areal_example_test.go b/stat/spatial_areal_example_test.go
new file mode 100644
index 0000000000..b1222da8b2
--- /dev/null
+++ b/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
+}
diff --git a/stat/spatial_example_test.go b/stat/spatial_example_test.go
new file mode 100644
index 0000000000..105bd8a0f2
--- /dev/null
+++ b/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
+}
diff --git a/stat/spatial_test.go b/stat/spatial_test.go
new file mode 100644
index 0000000000..b413cb75e1
--- /dev/null
+++ b/stat/spatial_test.go
@@ -0,0 +1,191 @@
+// 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"
+	"math/rand"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/matrix/mat64"
+)
+
+func simpleAdjacency(n, wide int, diag bool) *mat64.Dense {
+	m := mat64.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) *mat64.Dense
+
+	wantMoranI float64
+	wantZ      float64
+	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 TestMoranI(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)
+
+		m := NewMoran(data, locality)
+
+		gotI := m.I()
+		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)
+		}
+		gotZ := m.Z()
+		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)
+		}
+	}
+}
+
+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)
+
+		g := NewGetisOrd(data, locality)
+
+		if g.Len() != test.n {
+			t.Errorf("unexpected length: got:%d want:%d", g.Len(), test.n)
+			continue
+		}
+
+		const thresh = 6
+		var nseg int
+		segstart := -1
+		for i := 0; i < g.Len(); i++ {
+			gi := g.Gstar(i)
+			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 && g.Len()-segstart >= 5 {
+			nseg++
+		}
+		if nseg != test.wantSegs {
+			t.Errorf("unexpected number of significant segments for test %d: got:%d want:%d", ti, nseg, test.wantSegs)
+		}
+	}
+}