github.com/gopherd/gonum@v0.0.4/stat/spatial/spatial.go (about)

     1  // Copyright ©2017 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package spatial
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/gopherd/gonum/mat"
    11  	"github.com/gopherd/gonum/stat"
    12  )
    13  
    14  // TODO(kortschak): Implement weighted routines.
    15  
    16  // GetisOrdGStar returns the Local Getis-Ord G*i statistic for element of the
    17  // weighted data using the provided locality matrix. The returned value is a z-score.
    18  //
    19  //  G^*_i = num_i / den_i
    20  //
    21  //  num_i = \sum_j (w_{ij} x_j) - \bar X \sum_j w_{ij}
    22  //  den_i = S \sqrt(((n \sum_j w_{ij}^2 - (\sum_j w_{ij})^2))/(n - 1))
    23  //  \bar X = (\sum_j x_j) / n
    24  //  S = \sqrt((\sum_j x_j^2)/n - (\bar X)^2)
    25  //
    26  // GetisOrdGStar will panic if locality is not a square matrix with dimensions the
    27  // same as the length of data or if i is not a valid index into data.
    28  //
    29  // See doi.org/10.1111%2Fj.1538-4632.1995.tb00912.x.
    30  //
    31  // Weighted Getis-Ord G*i is not currently implemented and GetisOrdGStar will
    32  // panic if weights is not nil.
    33  func GetisOrdGStar(i int, data, weights []float64, locality mat.Matrix) float64 {
    34  	if weights != nil {
    35  		panic("spatial: weighted data not yet implemented")
    36  	}
    37  	r, c := locality.Dims()
    38  	if r != len(data) || c != len(data) {
    39  		panic("spatial: data length mismatch")
    40  	}
    41  
    42  	n := float64(len(data))
    43  	mean, std := stat.MeanStdDev(data, weights)
    44  	var dwd, dww, sw float64
    45  	if doer, ok := locality.(mat.RowNonZeroDoer); ok {
    46  		doer.DoRowNonZero(i, func(_, j int, w float64) {
    47  			sw += w
    48  			dwd += w * data[j]
    49  			dww += w * w
    50  		})
    51  	} else {
    52  		for j, v := range data {
    53  			w := locality.At(i, j)
    54  			sw += w
    55  			dwd += w * v
    56  			dww += w * w
    57  		}
    58  	}
    59  	s := std * math.Sqrt((n-1)/n)
    60  
    61  	return (dwd - mean*sw) / (s * math.Sqrt((n*dww-sw*sw)/(n-1)))
    62  }
    63  
    64  // GlobalMoransI performs Global Moran's I calculation of spatial autocorrelation
    65  // for the given data using the provided locality matrix. GlobalMoransI returns
    66  // Moran's I, Var(I) and the z-score associated with those values.
    67  // GlobalMoransI will panic if locality is not a square matrix with dimensions the
    68  // same as the length of data.
    69  //
    70  // See https://doi.org/10.1111%2Fj.1538-4632.2007.00708.x.
    71  //
    72  // Weighted Global Moran's I is not currently implemented and GlobalMoransI will
    73  // panic if weights is not nil.
    74  func GlobalMoransI(data, weights []float64, locality mat.Matrix) (i, v, z float64) {
    75  	if weights != nil {
    76  		panic("spatial: weighted data not yet implemented")
    77  	}
    78  	if r, c := locality.Dims(); r != len(data) || c != len(data) {
    79  		panic("spatial: data length mismatch")
    80  	}
    81  	mean := stat.Mean(data, nil)
    82  
    83  	doer, isDoer := locality.(mat.RowNonZeroDoer)
    84  
    85  	// Calculate Moran's I for the data.
    86  	var num, den, sum float64
    87  	for i, xi := range data {
    88  		zi := xi - mean
    89  		den += zi * zi
    90  		if isDoer {
    91  			doer.DoRowNonZero(i, func(_, j int, w float64) {
    92  				sum += w
    93  				zj := data[j] - mean
    94  				num += w * zi * zj
    95  			})
    96  		} else {
    97  			for j, xj := range data {
    98  				w := locality.At(i, j)
    99  				sum += w
   100  				zj := xj - mean
   101  				num += w * zi * zj
   102  			}
   103  		}
   104  	}
   105  	i = (float64(len(data)) / sum) * (num / den)
   106  
   107  	// Calculate Moran's E(I) for the data.
   108  	e := -1 / float64(len(data)-1)
   109  
   110  	// Calculate Moran's Var(I) for the data.
   111  	//  http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-how-spatial-autocorrelation-moran-s-i-spatial-st.htm
   112  	//  http://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-global-morans-i-additional-math.htm
   113  	var s0, s1, s2 float64
   114  	var var2, var4 float64
   115  	for i, v := range data {
   116  		v -= mean
   117  		v *= v
   118  		var2 += v
   119  		var4 += v * v
   120  
   121  		var p2 float64
   122  		if isDoer {
   123  			doer.DoRowNonZero(i, func(i, j int, wij float64) {
   124  				wji := locality.At(j, i)
   125  
   126  				s0 += wij
   127  
   128  				v := wij + wji
   129  				s1 += v * v
   130  
   131  				p2 += v
   132  			})
   133  		} else {
   134  			for j := range data {
   135  				wij := locality.At(i, j)
   136  				wji := locality.At(j, i)
   137  
   138  				s0 += wij
   139  
   140  				v := wij + wji
   141  				s1 += v * v
   142  
   143  				p2 += v
   144  			}
   145  		}
   146  		s2 += p2 * p2
   147  	}
   148  	s1 *= 0.5
   149  
   150  	n := float64(len(data))
   151  	a := n * ((n*n-3*n+3)*s1 - n*s2 + 3*s0*s0)
   152  	c := (n - 1) * (n - 2) * (n - 3) * s0 * s0
   153  	d := var4 / (var2 * var2)
   154  	b := d * ((n*n-n)*s1 - 2*n*s2 + 6*s0*s0)
   155  
   156  	v = (a-b)/c - e*e
   157  
   158  	// Calculate z-score associated with Moran's I for the data.
   159  	z = (i - e) / math.Sqrt(v)
   160  
   161  	return i, v, z
   162  }