Implementing Neatest Neighborhood Gaussian Process (NNGP) in Stan

[Lu Zhang]

Hi Stan,

      I am working on the implementation of Neatest Neighborhood Gaussian Process (NNGP) in Stan. I wrote a test folder and here is the share link: https://drive.google.com/drive/folders/0ByIcksYvjph9aS1GajdzY2dCWGc?usp=sharing.

The folder exceeds the limits for this group so I put it in the google driver. Please let me know if it is not accessible. It contains the synthetic data and codes for four different models. For more information, please check "readme" and "summary.pdf". Since it hasn't been completed yet, the summary report is still in editing. Welcome for suggestions for improving my code. I also posted two questions in the end of this topic. Thank you very much!

The following is file "Readme.txt" which gives a brief introduction of this folder:

# Introduction:

This is a test folder for implementing Full Gaussian process and Nearest 

Neighborhood Gaussian Process on spatial data in Stan. For more details 

about synthetic data, models, and simulation results, please refer to 

“summary.pdf”.  

# How to run the program:

The language interface to perform the analysis is Rstan. All the executable files 

are in subfolder “projects” and Stan codes are saved in subfolder “src”. To run 

the analysis, open R program in “projects” and set the working directory to be 

the path of this folder. R code in subfolder “projects” will load data and link 

Stan codes with a relative path. 

# Folder Structure:

data: Synthetic data and parameters for priors

projects:         Rstan codes for two models.

results:         Simulation results

src:         Stan codes

The following are two questions I met, more information can be found in "summary.pdf":

1. In model1, the performance of FullGP highly depends on initial values. NNGP did not have such a high number of divergent transitions with the same initials and number of iterations as follows. What's more, it converged. I want to compare FullGP and NNGP with initial values far from the true value, but the results from FullGP didn't even converge.  

Initial value: beta = c(0, 0), sigmasq = 6, tausq = 1, phi= 10;

With 400 iterations, I got Warning messages:

1: There were 199 divergent transitions after warm up. Increasing adapt_delta above 0.8 may help. 

2: Examine the pairs() plot to diagnose sampling problems

2. The Full Gaussian Process in model 2 takes a very long time to run. I got 1200 out of 2000 iterations over 10 hours and then my computer stopped running. Hope someone could give me some suggestions on the coding.

Best Wishes,

Lu Zhang

[Mike Lawrence]

Regarding your question #2, on improving the efficiency of the fullGP models:

First, in your custom function "get_w_L" you should take advantage of the fact that the covariance matrix is symmetric by doing:

    for (j in 1:(N-1)){

        for (k in j:N){

            temp_dist[j, k] = sigmasq * exp(- phi * dist[j, k]) ;

            temp_dist[k, j] = temp_dist[j, k] ;

        }

    }

  for(i in 1:N){

  temp_dist[i,i] = sigmasq + tausq ;

  }

You might also be getting into trouble with your priors; for GPs, I follow the guidelines posted at the bottom of the wiki: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations

--
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.

[Rob Trangucci]

This might be a little faster (for any function that builds the GP covariance matrix):

for (j in 1:(N-1)) {

       temp_dist[j, j] = sigmasq + tausq;

        for (k in (j+1):N) {

            temp_dist[j, k] = sigmasq * exp(- phi * dist[j, k]);

            temp_dist[k, j] = temp_dist[j, k];

        }

    }

temp_dist[N,N] = sigmasq + tausq;

Rob

[Rob Trangucci]

And for model fullGP_w.stan, I would rewrite the parameters/model block like follows:

...

parameters {

    vector[2] beta;

    real<lower = 0> sigmasq;

    real<lower = 0> tausq;

    real<lower = ap, upper = bp> phi;

    vector[N] eta;

}

model {

    vector[N] w;

    matrix[N, N] L_cov;

    L_cov = get_w_L(sigmasq, tausq, phi, dist, N);

    sigmasq ~ inv_gamma(as, bs);

    tausq ~ inv_gamma(at, bt);

    eta ~ normal(0, 1);

    w = mu_w + L_cov * eta;

    for (i in 1:N){

        Y[i] ~ normal((X[i,] * beta + w[i]), tausq);

    }

}

Rob

[Mike Lawrence]

Oh, and just noticed you have flat priors on your covariates; have you tried anything more informative? Probably not the cause of the difference in performance between FullGP & NNGP, but making them even a little informative should help reduce the compute time for both (let alone solid arguments for non-sensicality of flat priors made previously on this list).

[Andrew Gelman]

Also I'd avoid that inverse-gamma prior, and I'd parameterize in terms of sigma and tau, not sigmasq and tausq.

A

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

[Lu Zhang]

I see, Thank you very much! I will change my code accordingly.

Lu

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

> On Sep 10, 2016, at 8:01 PM, Rob Trangucci <robert.t...@gmail.com> wrote:
>
> And for model fullGP_w.stan, I would rewrite the parameters/model block like follows:
> ...
> parameters {
> vector[2] beta;
> real<lower = 0> sigmasq;
> real<lower = 0> tausq;
> real<lower = ap, upper = bp> phi;
> vector[N] eta;
> }
> model {
> vector[N] w;
> matrix[N, N] L_cov;
> L_cov = get_w_L(sigmasq, tausq, phi, dist, N);
> sigmasq ~ inv_gamma(as, bs);
> tausq ~ inv_gamma(at, bt);
> eta ~ normal(0, 1);
> w = mu_w + L_cov * eta;
> for (i in 1:N){
> Y[i] ~ normal((X[i,] * beta + w[i]), tausq);
> }
> }

You could vectorize the loop to

Y ~ normal(X * beta + w, tau);

As long as Y is a size-N vector and X is an N x 2 matrix.

And rather than defining tau = sqrt(tausq), you can
just make tau a parameter and put a prior on it
directly.

- Bob

[Bob Carpenter]

Also, before running 2000 iterations, try running 200 and
see what you get. Then 400, and so on, until the convergence
diagnostic R-hat is close to 1 and the effective sample size
you need.

- Bob

> 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.
>
>
> --
> 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.
>
>
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.

[Lu Zhang]

Thank you for your advice. Yes, I tried very informative priors for regression coefficients in the very beginning and Stan gave me good results. I may try a multivariate normal with large variance later. But aside from the issue of convergence, will flat prior increase computing time?

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Lu Zhang]

Thank you! Actually, the parameter I put is tausq, not tau. I am not very clear why I should set a prior for tau instead of tausq. Is there any beneficial I can get from doing that?

 

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Lu Zhang]

I see, thank you very much! 

> 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.

> 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 a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

> To unsubscribe from this group and all its topics, 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.

> 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 a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Andrew Gelman]

Which is something that at some point we'd like to automate.

[Andrew Gelman]

1.  Inv-gamma has problems (see Gelman, 2006).

2.  tau is more interpretable than tausq.  You can have prior info on tau, it's not so clear how to do that with tausq

3.  Normal dist in Stan uses the normal(mean, sd) parameterization, NOT normal(mean, var).

That pretty much covers it!

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.

[Lu Zhang]

I see! Thank you so much!

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.

[Lu Zhang]

Now I see a lot of problems in my codes, I will modify my code tonight and update my code later. Thank you very much!

Best Wishes,

Lu Zhang

>> 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.

>> 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 a topic in the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

>> To unsubscribe from this group and all its topics, 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.

>> 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 a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Yes flat prior will probably increase computing time
relative to a weakly informative or informative prior.

The advantage of putting a prior on the scale (std dev) rather
than squared scale (variance) is that the scale parameter has the
same units as the mean, making it easier and more natural to work
with. A perhaps bigger issue is the behavior of the prior as the
scale approaches zero and how much mass it spreads out into
the tail.

- Bob

> 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.
>
>
> --
> 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.
>
>
>
> --
> 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.
>
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.

[Lu Zhang]

Hi Bob, 

    Thanks, I tested a multivariate normal prior for intercept and slope. It converges faster and solves the first problem I posted in the topic. As for priors for variances, I set a normal prior in FullGP for tau, but I still use the inverse-Gamma for "sigmasq" so that I don't need to take sigma^2 later. What's more, I think variogram is more widely used in spatial analysis, the mean of the prior for "sigmasq" is actually got from the empirical variogram in the exploration of data analysis ( please check EDA in subfolder "project"). So I think it is also nature to work with "sigmasq."

In fullGP.stan, I changed my code into:

...

...

parameters {

    vector[2] beta;

    real<lower = 0> sigmasq;

    real<lower = 0> tausq;

    real<lower = ap, upper = bp> phi;

}

model{

    beta ~ multi_normal(uB, VB);

    sigmasq ~ inv_gamma(as, bs);

    tausq ~ inv_gamma(at, bt);

    Y ~ multi_normal_cholesky(X * beta, fullGP_L(sigmasq, tausq, phi, dist, N));

}

So should I take the Cholesky decomposition of the covariance matrix (VB) in the prior of beta in "transformed data", and use multi_normal_cholesky() instead of multi_normal()? 

Thanks,

Lu Zhang

> 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.

> 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 a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

> To unsubscribe from this group and all its topics, 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.

> 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 a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Lu Zhang]

Hi Bob,

    I have a question about coding in FullGP_w.stan. In your previous email, you suggested a rewrite of my parameter/model block for FullGP_w.stan: 

> ...
> parameters {
>     vector[2] beta;
>     real<lower = 0> sigmasq;
>     real<lower = 0> tausq;
>     real<lower = ap, upper = bp> phi;
>     vector[N] eta;
> }
> model {
>     vector[N] w;
>     matrix[N, N] L_cov;
>     L_cov = get_w_L(sigmasq, tausq, phi, dist, N);
>     sigmasq ~ inv_gamma(as, bs);
>     tausq ~ inv_gamma(at, bt);
>     eta ~ normal(0, 1);
>     w = mu_w + L_cov * eta;
>     for (i in 1:N){
>         Y[i] ~ normal((X[i,] * beta + w[i]), tausq);
>     }
> }

Different from FullGP.stan, FullGP_w.stan includes the spatial random effect "w" in parameter space, so the dimension increases dramatically from 5 to 1005.   But if I use the code above, how should I retrieve "w"? Should I use a transformed parameter block for that? Thank you very much!

Best Wishes,

Lu  

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

I'm not the right person to answer in-depth GP questions.

- Bob

> To unsubscribe from this group and all its topics, 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.
>
>

> --
> 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.

[Seth Flaxman]

Hi Lu. Yes, use multi_normal_cholesky() instead of multi_normal()
wherever possible. And yes, if you put w in the transformed_parameter
block it'll be available for you to use other places.

Seth

[Lu Zhang]

Hi Seth. Thank you very much! I will change accordingly.

>> To unsubscribe from this group and all its topics, 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.

>> 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 a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JwmsztFcqWM/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Lu Zhang]

Hi Stan,

    I just updated my test folder. Here are the changes in comparison to the previous version:

# Updates:

1. Input lower triangular part of the distance matrix instead of the whole distance matrix.

2. Change the priors of Beta from flat into multivariate normal.

3. Change the initial value into 

    Model 1:  beta (0,0); sigmasq = 4; tausq = 1; phi = 10;  

    Model 2:  beta (0,0); sigma = 2; tau = 1; phi = 10; 

4. Change the priors of sigma and tau in model 2 into half-Cauchy

5. Use a centred parametrisation for beta1 and w instead of an uncentered w in model 2 

6. Add all the simulation results in subfolder “results”

7. Add summarization of all the simulations in “summary.pdf”

Thanks,

[Craig Wang]

Hi Lu,

First of all, thank you for implementing and sharing your code with good documentation. I would like to ask you, why do you switch to a centred (at beta1) parameterization for w, the Stan manual mentioned non-centred is preferred?

[Lu Zhang]

Hi Craig,

    Sorry for the late reply. I use a centred parameterization because I observed a high correlation between my intercept and my zero centred random effect. You can look a paper in Biometrika by Gelfand, Sahu and Carlin about hierarchical centering. Hope this helps.

Best Wishes,

Lu Zhang

[Craig Wang]

Hi Lu, 

Thanks for directing me to the paper, I will have a read.

Best,

Craig

To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.