# Estimating parameters for a spatial process

I'm given an $n\times n$ grid of positive integer values. These numbers represent an intensity that should correspond to the strength of belief of a person occupying that grid location (a higher value indicating a higher belief). A person will in general have an influence over multiple grid cells.

I believe that the pattern of intensities should "look Gaussian" in that there will be a central location of high intensity, and then the intensities taper off radially in all directions. Specifically, I'd like to model the values as coming from a "scaled Gaussian" with a parameter for the variance and another for the scale factor.

There are two complicating factors:

• the absence of a person will not correspond to a zero value, because of background noise and other effects, but the values should be smaller. They can be erratic though, and at a first approximation might be difficult to model as simple Gaussian noise.
• The intensity range can vary. For one instance, the values might range between 1 and 10, and in another, between 1 and 100.

I'm looking for an appropriate parameter estimation strategy, or pointers to relevant literature. Pointers to why I'm approaching this problem the wrong way altogether would also be appreciated :). I've been reading about kriging, and Gaussian processes, but that seems like very heavy machinery for my problem.

• What do you mean by a Gaussian with a variance and scale parameter? The variance parameter is the scale parameter of a Gaussian! I'm also a little unsure about the model you've set up so far. Can you describe the problem you're actually trying to solve in more detail? Using a Gaussian to model low granularity integer-valued observations seems fishy. – cardinal Apr 7 '11 at 12:19
• (+1) For an interesting question. Looking forward to understanding what you're trying to solve a bit better. – cardinal Apr 7 '11 at 12:22
• Here are several observations: 1. If your values are integer, using Gaussian does not seem appropriate. 2. It is not clear what is the purpose of your model, do you want to identify the clusters of strong belief for example? What will be interpretation of your parameters if you had them? 3. Since you have a grid, why not try to fit mixture of bivariate distributions? Then the grid will be the support of the distribution (say unit square) and the intensities will correspond to high probability regions. – mpiktas Apr 7 '11 at 12:28
• Thanks all for the interesting points. Let me attempt to clarify. The choice of "Gaussian", in the light of comments, might be a red herring that causes more confusion than it helps. The key feature of the data is the high intensity values at point of highest belief in the location of the person, and the tapering off "radially" around it (which I've observed empirically). The intensity values come from the solution to a (linear) inverse problem, and so actually don't necessarily need to be integral - that's just the data we happen to have. – Suresh Venkatasubramanian Apr 7 '11 at 22:00
• btw I appreciate the attempts to make the question more well defined and better modelled. I'll do my best to explain the actual data setting so as to converge on the right modelling assumptions. – Suresh Venkatasubramanian Apr 7 '11 at 22:03

You can use this module of the pysal python library for the spatial data analysis methods I discuss below.

Your description of how each person's attitude is influenced by the attitudes of the people surrounding her can be represented by a spatial autoregressive model (SAR) (also see my simple SAR explanation from this SE answer 2). The simplest approach is to ignore other factors, and estimate the strength of the influence how surrounding people affect one another's attitudes by using the Moran's I statistic.

If you want to assess the importance of other factors while estimating the strength of the influence of surrounding people, a more complex task, then you can estimate the parameters of a regression: $y = bx + rhoWy + e$. See the docs here.(Methods of estimating this type of regression come from the field of spatial econometrics and can get much more sophisticated than the reference I gave.)

Your challenge will be to build a spatial weights matrix ($W$). I think each element $w_{ij}$ of the matrix should be 1 or 0 based on whether the person $i$ is within some distance you feel that it is required to influence the other person $j$.

To get an intuitive idea of the problem, below I illustrate how a spatial autoregressive data generating process (DGP) will make a pattern of values. For the 2 lattices of simulated values the white blocks represent high values and the dark blocks represent low values.

In the first lattice below the grid values have been generated by a normally distributed random process (or Gaussian), where $rho$ is zero. In the next lattice below the grid values have been generated by a spatial autoregressive process, where $rho$ has been set to something high, say .8. • That's very interesting (and so is the related Geary C). This might be close to what I need. – Suresh Venkatasubramanian May 10 '11 at 19:16
• Geary C helps you see how values close to one another cluster, even values in the middle of the distribution. Moran's I helps you see how very high values cluster with very high values, and very low values cluster around very low values. So perhaps you are correct and the simplest and best method is Geary's C. Remember the Geary's C approach is exploratory and will not let you condition your results on other factors. Look at this python module for code to run Geary's C: pysal.org/1.1/library/esda/geary.html. – b_dev May 10 '11 at 20:57
• Let me play with these some more. If it seems to do what I need (and I think it will), this sounds like the best answer. – Suresh Venkatasubramanian May 10 '11 at 21:13

Here is a simple idea which might work. As I've said in the comments if you have a grid with intensities why not fit density of bivariate distribution?

Here is the sample graph to illustrate my point: Each grid point with is displayed as a square, colored according to intensity. Superimposed on the graph is the contour plot of bivariate normal density plot. As you can see the contour lines expand in the direction of decreasing intensity. The center will be controled by the mean of bivariate normal and the spread of the intensity according to covariance matrix.

To get the estimates of mean and covariance matrix simple numerical optimisation can be used, compare the intensities with values of density function using the mean and the covariance matrix as parameters. Minimize to get the estimates.

This is of course strictly speaking not a statistical estimate, but at least it will give you an idea how to proceed further.

Here is the code for reproducing the graph:

require(mvtnorm)
sigma=cbind(c(0.1,0.7*0.1),c(0.7*0.1,0.1))

x<-seq(0,1,by=0.01)
y<-seq(0,1,by=0.01)
z<-outer(x,y,function(x,y)dmvnorm(cbind(x,y),mean=mean,sigma=sigma))

mz<-melt(z)

mz$X1<-(mz$X1-1)/100
mz$X2<-(mz$X2-1)/100

colnames(mz)<-c("x","y","z")

mz$intensity<-round(mz$z*1000)

ggplot(mz, aes(x,y)) + geom_tile(aes(fill = intensity), colour = "white") + scale_fill_gradient(low = "white",     high = "steelblue")+geom_contour(aes(z=z),colour="black")


Your model is a two-dimensional random field $X[i,j]$, and you are trying to estimate the joint distribution of the integer-values random variables $X[i,j]$. You will want to assume spatial stationarity: that is, the joint distribution of $(X[i_1,j_1],...,X[i_m,j_m])$ is the same as the joint distribution of $(X[i_1+k,j_1+l]...,X[i_m+k,j_m+l])$. In particular, the marginal distribution is the same for every cell. A simple question to ask is the autocorrelation structure of the field. That is, what is $corr(X[i_1,j_1],X[i_2,j_2])$ given the distance $d([i_1,j_1],[i_2,j_2])$? We represent this as a function $\rho(d)$. A simple model for the autocorrelation structure is $\rho(d)=kd^{-1}$, where $k$ is a constant.

A 'gaussian' effect corresponds to a quadratic distance function, but there are many other distance functions you should consider, such as the taxicab norm $d([i_1,j_1],[i_2,j_2]) = |i_1-i_2|+|j_1-j_2|$. Once you have decided on a distance function and the form of your model for autocorrelation it is simple enough to estimate $\rho(d)$ e.g. via maximum likelihood. For more ideas, look for "random field".

• "Want to assume spatial stationarity" appears directly to contradict the OP's assumption that "the intensities taper off radially in all directions." – whuber Apr 7 '11 at 19:52
• How so? Such a pattern would occur with the autocorrelation structure I proposed. – charles.y.zheng Apr 7 '11 at 19:57
• @charles It's an important point: if indeed this apparent trend is to be attributed to autocorrelation, then in principle another independent realization of the process could appear to have a dramatically different trend, such as an increase in value away from a central point. Because the OP has clearly articulated and distinguished some deterministic elements to the trend (the "radial tapering") and correlational elements ("have influence over multiple grid cells"), a reply that respects this would likely be viewed more positively than one that asserts the OP "will want" to change his mind. – whuber Apr 7 '11 at 20:04
• I'm not sure I understand the spatial stationarity condition. On the surface, it seems to be at odds with the idea of having a "peak that tapers out" at a specific location, but I'm clearly not understanding something. – Suresh Venkatasubramanian Apr 7 '11 at 22:02
• @charles, the pattern you describe will be present for each of grid points, due to spatial stationarity assumption. Stationarity is basically saying that all my points behave similarly. This is not the case described by OP. The answer still is very good, but not appropriate in this case. – mpiktas Apr 8 '11 at 11:28