1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
// 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"
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/stat"
)
// TODO(kortschak): Implement weighted routines.
// GetisOrdGStar returns the Local Getis-Ord G*i statistic for element 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
if doer, ok := locality.(mat.RowNonZeroDoer); ok {
doer.DoRowNonZero(i, func(_, j int, w float64) {
sw += w
dwd += w * data[j]
dww += w * w
})
} else {
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)
doer, isDoer := locality.(mat.RowNonZeroDoer)
// Calculate Moran's I for the data.
var num, den, sum float64
for i, xi := range data {
zi := xi - mean
den += zi * zi
if isDoer {
doer.DoRowNonZero(i, func(_, j int, w float64) {
sum += w
zj := data[j] - mean
num += w * zi * zj
})
} else {
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
if isDoer {
doer.DoRowNonZero(i, func(i, j int, wij float64) {
wji := locality.At(j, i)
s0 += wij
v := wij + wji
s1 += v * v
p2 += v
})
} else {
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
}
|
