Priors for spatial lengthscale

[Christopher Peterson]

I'm trying to model the mean values of data collected from a number of sites where means at each site can be spatially autocorrelated.  I used a squared-exponential covariance function to build the spatial correlation matrix.Unfortunately, the model has a large number of diverging errors; it also fails with the cholesky decomposition of the matrix frequently enough for me to think that it's a substantial problem.  (Message: "Exception thrown at line 36: cholesky_decompose: Cholesky decomposition of m failed")

This model has no difficulties if the spatial component is removed. I think the most likely source of the error would be the prior for the spatial lengthscale.  I'd like for this parameter to be able to shrink to near-0 (which result in an identity matrix for the spatial correlation) while still having sufficient mass in the tail to allow for substantial spatial signal.  I've tried a variety of options (half Cauchy, half student_t, exponential), but none seem to have improved anything. 

Is there something else wrong with this model?  Should I be using a different covariance function? Does anyone have any suggestions or advice for the prior?

Thank you

/* Model Objectives:
Estimate the mean value at each site and the grand mean.
Allow for standard deviation to vary among sites

Account for spatial autocorrelation of means among sites
*/

data{

  int n; // observations

  int nSite; // sites

  int siteID[n]; // which site belongs to each obs. 

  vector[n] x;

  matrix[nSite,nSite] Dist; // Distance matrix, pre-scaled so the minimum distance is 1

}

transformed data{

  real xScale; // the range of x divided by 4, to establish a reasonable scale for priors

  real xMean;

  matrix[nSite,nSite] sqDist; // pre-compute squared distance matrix

  xMean <- mean(x);

  xScale <- (max(x) - min(x))/4;

  for(i in 1:nSite) for(j in 1:nSite)

    sqDist[i,j] <- Dist[i,j]^2;

}

parameters{

  real<lower=0> spatialScale; // spatial lengthscale; 

  real grandMean; // overall mean of x

  

  // Non-central parameterizations; use is defined in next block

  vector[nSite] raw_siteMean;

  real<lower=0,upper=pi()/2> raw_siteSD;

  real<lower=0,upper=pi()/2> raw_indSDHier;

  vector<lower=0,upper=pi()/2>[nSite] raw_indSD;

}

transformed parameters{

  vector[nSite] siteMean; // mean of x at each site

  real siteSD; // stdev of siteMean

  real indSDHier; // overall stdev of x around the sitemeans

  vector[nSite] indSD; // each site's scaling of indSDHier

  

  matrix[nSite, nSite] Sigma_spatial; // Spatial matrix

  matrix[nSite,nSite] L_Sigma_spatial; // cholesky of spatial matrix

  Sigma_spatial <- exp(-sqDist/(2*spatialScale^2)); // squared exponential covariance

  L_Sigma_spatial <- cholesky_decompose(Sigma_spatial);

  

  siteSD <- tan(raw_siteSD) * xScale; // ~half-cauchy(0, xScale)

  siteMean <-  siteSD * 

    L_Sigma_spatial * // comment out this line to remove spatial dependence

    raw_siteMean + grandMean; // ~ MVN_chol(grandMean, siteSD*L_Sigma_spatial); 

  indSDHier <- tan(raw_indSDHier) * xScale; // ~half-cauchy(0, xScale)

  for(i in 1:nSite)

    indSD[i] <- indSDHier * tan(raw_indSD[i]); // ~half-cauchy(0, indSDHier)

}

model{

  grandMean ~ normal(xMean,xScale);

  raw_siteMean ~ normal(0,1);

  spatialScale  ~ cauchy(0, 1); // A variety of other priors have been attempted, including student_t and exponential

  x ~ normal(siteMean[siteID], indSD[siteID]);

}

 

If it helps, I can provide some data and a set-up script. This was all run in rstan 2.9.

[Michael Betancourt]

The length scale shrinking to zero does not give

you an identity covariance but rather a covariance

where all entries are unity, which is not positive-definite

and hence will frustrate the Cholesky decomposition.

A prior will still be helpful, but you really need to add

a multiple of the identity matrix, with the scaling another

possible parameter, to Sigma_spatial to regularize the

covariance.

On Jun 21, 2016, at 9:07 PM, Christopher Peterson <christopher...@gmail.com> wrote:

I'm trying to model the mean values of data collected from a number of sites where means at each site can be spatially autocorrelated.  I used a squared-exponential covariance function to build the spatial correlation matrix.Unfortunately, the model has a large number of diverging errors; it also fails with the cholesky decomposition of the matrix frequently enough for me to think that it's a substantial problem.  

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mike Lawrence]

Christopher, have you tried using an inverse-lengthscale parameterization? You simply change this:

    Sigma_spatial <- exp(-sqDist/(2*spatialScale^2)); 

to

    Sigma_spatial <- exp(-sqDist*(2*spatialScale^2)); 

I find this a more intuitive parameterization as values from 0-1 imply near-monotonicity within the range of the data (if x is scaled to 0-1), making the typical peaked-at-zero priors sensibly reflecting a bias for simpler functions. 

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[a...@ariddell.org]

It would be wonderful to capture some of this discussion and put it on
the Stan Wiki, somewhere on
https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations .
Right now the page does not have any discussion of priors for GPs.

Best,

A

[Avraham Adler]

Honestly, there are so many gems sprinkled through these hundreds of threads, that it may pay for someone to collate and archive them.

Anyone interested in editing the inaugural issue of "Proceedings of the General Bayesian Discussions & Friends of Stan"? 8-)

Avi

> > email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> >
> >
> > --
> > You received this message because you are subscribed to the Google Groups
> > "Stan users mailing list" group.
> > To unsubscribe from this group and stop receiving emails from it, send an

> > email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> >
>
> --
> You received this message because you are subscribed to the Google Groups
> "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an

> email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

Honestly we’re just throwing this stuff out as best we can

given the inspiration of the original question and we (the

dev team and other common answerers) really are the 

worst to judge what’s useful for users of various background

and what’s not.  But it would be great if users would collate

what they found useful to make it easier for others to find

and give us feedback as to what explainations/recommendations

worked well.  And it would be doubly great if someone built

some site or forum to facilitate users identifying/sharing 

anything they find useful.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I agree that we should keep building out the Wiki.

The real problem we have is finding things. Google groups
seem to be a very bad forum for that compared to Stack Overflow.
But then some of our questions are more of the Cross Validated
variety.

- Bob

[Mike Lawrence]

Done! :O)

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Mike Lawrence]

Bob, have the devs considered discourse as a mailing-list format replacement? https://www.discourse.org/about/

[Daniel Lee]

1. Mind creating a new thread for separate discussion topics? One of the difficulties in finding things is due to threads getting hijacked.

2. Yes, we've considered discourse and flarum. There just hasn't been enough expertise / benefit / time to move from what sort of works to another system that has different flaws. If anyone wanted to step up, set it up, and help through the migration process, I think we'd love the help. Unfortunately, discourse fell through due to: cost and lack of some key features (that I can't remember now). But really, cost. In both time and $.

[Andrew Gelman]

I suppose we could pay someone to archive these discussions, but if in the meantime one of you could add something to the wiki on priors, that would be great!

As Alan said, it's' here:

https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations 

Andrew

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Mike Lawrence]

I added some recommendations regarding GP parameters, including the one I made in this thread, to the bottom of the wiki earlier today. Given my own novice standing with regards to GPs, it might be beneficial for someone to check my reasoning/understanding.

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Christopher Peterson]

Regarding the original question:

I'm a bit confused as to how that would work.  Given the squared-exponential covariance function Cov[i,j] = sigma^2 * exp(-D[i,j]^2 / l^2), where D[i,j] is the fixed distance between points i and j and l is the length-scale parameter, I was under the impression that as l^2 approached 0,  exp(-D[i,j]^2 / l^2) would also approach 0 for the off-diagonals (while the diagonals, with 0 distance, remained at 1).  Thus, a prior that shrunk l towards 0 should reduce the correlation to an identity matrix.  Is this mistaken?  I know I must be screwing something up somewhere, but I don't see how this is it. 

Thank you for your help, everyone.

[Andrew Gelman]

Thanks. 

[Mike Lawrence]

Here's some R code to help intuition:

  library(MASS)

  inv_rho = .001 #contrast .001 & 10

  eta = 1

  sigma = 1

  N = 1e2

  s = matrix(NA,nrow=N,ncol=N)

  for(i in 1:N){

    for(j in 1:N){

      s[i, j] = eta * exp(-((i-j)^2)*inv_rho) + sigma

    }

  }

  out = MASS::mvrnorm(1,rep(0,times=N),s)

  plot(out,type='l',ylim=c(-10,10))

As has been pointed out in other threads here, inv_rho & eta play somewhat similar roles in areas of the parameter space, but generally you should observe that for a given value of eta, you get more wiggly plots with larger values for inv_rho.

[Andrew Gelman]

I think that informative priors on lenghscale can make a lot of sense.  To start with:  at the extremes of lengthscales near 0 or infinity, there will be no info to learn about this hyperparameter.  And, a lot of the time, you can get some (soft) bounds on reasonable lenghtscale based on the nature of the problem.  For example if we're modeling annual economic data, we probably wouldn't want a lenghtscale less than 2 years.  And so on.

A

[Mike Lawrence]

My bad, the code should be:

  library(MASS)

  inv_rho = exp(5) #compare exp(0) to exp(10)

  eta = 1

  sigma = .001 #setting this small to highlight effect of inv_rho

  N = 1e3

  s = matrix(NA,nrow=N,ncol=N)

  x = seq(0,1,length.out=N)

  for(i in 1:N){

    for(j in 1:N){

      s[i, j] = eta * exp(-((x[i]-x[j])^2)*inv_rho)

    }

  }

  for(i in 1:N){

    s[i, i] = eta + sigma

  }

  y = MASS::mvrnorm(1,rep(0,times=N),s)

  plot(x,y,type='l',ylim=c(-10,10))

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~