gaussian process fitting and posterior predictive distribution in stan

[Linas Mockus]

Hi,

I am exploring the possibility of simultaneous fitting and posterior predictive distribution in stan and combined the code from sections 16.3 and 16.4 of stan user's guide.  The combined code is attached.  I wonder if there are any possibilities to improve the code in order to get predictive posterior distribution faster (such as Cholesky Factorization).  Perhaps I am not able to use the analytical form outlined in 16.4 because eventually I will use informative priors for sigma_sq.

Can you share any insights into this?

Thanks,
Linas

[Bob Carpenter]

You can certainly Cholesky factor. Then you have to do
your own transform, but I thought there was an example of that in
the manual.

You should probably also look about doing anything predictive in
the generated quantities block. There’s a multi_normal_rng, but not
yet a multi_normal_prec_rng or multi_normal_cholesky_rng.

- Bob

> --
> 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.
> For more options, visit https://groups.google.com/groups/opt_out.
> <fitpredpost.stan>

[Linas Mockus]

Bob,

Thanks for the suggestion. 

Unfortunately while the GP fit converged the GP fit+prediction started to diverge (without mentioning that it is slow...).  I wonder if it is possible to fit using stan and then do the prediction inside R using stan posteriors of GP hyperparameters.  Can anybody suggest an example?  I am still novice in GP/Bayesian landscape.  Attached are stan codes to do the fitting using GP (borrowed from stan user's manual) and the code to do both fit and prediction with corresponding trace and simulated data.

I appreciate any help,
Linas

[Bob Carpenter]

You certainly can do that, but it's not going to give
you a different answer than doing it within Stan (unless there's
a bug on one side or the other). Stan's extract function lets
you pull out the parameters for each draw in the sample.

I'm guessing from your traceplot, which shows a huge increase
in variance after adaptation, that what's happened is that the
last step of adaptation found a smaller step size and was thus
able to explore more of the posterior. That is, the post-warmup
draws showing higher variance may be closer to the posterior
than the warmup iterations. My suggestion would be to run
warmup for more iterations and see what the posterior looks like.

You want to inspect what happened to adaptation in the two chains you're
showing. There's no way the posterior measures of convergence R-hat
and n_eff will look good given that traceplot with two chains of very
different means.

- Bob

> --
> 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.
> For more options, visit https://groups.google.com/groups/opt_out.

> <FitPredict.jpg><FitPredict.stan><Fit.jpg><Fit.stan><SimulatedData.jpg>

[Linas Mockus]

Bob,

Thanks for suggestion.  I will certainly try it.

My question is more why the traceplot for fit code shows very good mixing while the traceplot for fit & predict code shows bad mixing.  One possibility is that stan code for fit & predict is not right since stan code from fit & predict is from stan user's manual.

Another thing I am looking is an example in R of GP prediction using parameter estimates extracted from stan.

Linas

[Bob Carpenter]

Any prediction going on in the generated quantities block
won't affect the actual sampling algorithm.

It's always good to have an independent test of code.

The thing I'd be most worried about is the difference in behavior
during adaptation and sampling. It's a sign that adaptation hadn't
fully adapted.

- Bob

[Linas Mockus]

Bob,

Yes, I agree that generated quantities block shouldn't affect the sampling. But the example in sec 16.4 of stan user's manual places data (y1) and posterior predictive (parameter y2) into one long vector inside model block.  I have modified example even further by trying to estimate eta, rho, and sigma by combining examples in 16.3 and 16.4.  I feel that I may be doing something wrong.  Therefore I also wanted independent verification using R but couldn't not find a decent example.

Linas 

[Bob Carpenter]

We don't have any examples that we've generated in R,
but it should be exactly the same thing. You'll just
need to bring in somehting like MASS to do multivariate
sampling and be careful about the difference between precision
and covariance matrix parameterizations.

- Bob

[Linas Mockus]

Thank you.  Does it mean that the code in stan combined by me for fitting & predicting is good and I just need a longer warmup?  If so I will try to replicate the code in R and post so other novices such as me could use it.

Linas

[Bob Carpenter]

I'm not sure that's the problem, but you should definitely use
a longer warmup period to make sure everything converges.

- Bob

[Linas Mockus]

Bob,

You were right, after increasing warmup chains started to mix.

In my excitement I have crafted stan code with generated quantities
block. The second code was much faster. However, the predictions are
off for the stan code with the generated quantities block (as you have
suggested previously). Both codes generate comparable posteriors. I
guess there is some subtle difference between both codes I am unable to
grasp. Can you please help?

Linas

[Bob Carpenter]

I'd suggest taking smaller steps in moving from your Stan1.stan
to Stan2.stan model. In particular, starting from Stan1.stan,
I'd suggest replacingn

for (n in 1 : N1) y [n] <- y1 [n];
for (n in 1 : N2) y [N1 + n] <- y2 [n];
y ~ multi_normal (mu, Sigma);

with just

y1 ~ multi_normal(mu,Sigma);

and verifying that mu and Sigma fit as intended, and then
moving the sampling for y2 to the generated quantities block.

- Bob

> --
> 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.
> For more options, visit https://groups.google.com/groups/opt_out.

> <Stan1.stan><Stan2.stan>

[andy...@cubist-asia.com]

Hi Linas/Bob,

This thread has been extremely helpful for me to get up to speed on fitting and predicting GP "in one go". Thanks very much guys!

Based on, FitPredict, I have tried to add multivariate input support and automatic relevance determination (Ch 14.2, 14.3, doc: 2.6.2)

1. As per suggestion in the manual, I have made the following changes:

parameters {
....
  vector<lower=0>[D] rho;
.....
}

model {
....
  matrix [N, N] Sigma;
  for (i in 1 : (N - 1)) {
    for (j in i : N) {
      vector[D] tmp_exp_term;
      tmp_exp_term <- rho .* (x[i] - x[j]);

      Sigma [i, j] <- eta_sq * exp (-dot_self(tmp_exp_term));   
    //    exp (-rho_sq * pow (x [i] - x [j], 2));
      Sigma [j, i] <- Sigma [i, j];
    }
  }
....
....
rho    ~ cauchy (0, 5);
}

I want to ask, from efficiency perspective, are there much better ways to do it without moving to Cholesky/analytical form?

2. Attached is a traceplot of selected paramters i have run, with N = 20, D=3. While the predicted y2 have converged quite well, trace of eta_sq and rho seems to have done less well. Is this something I should worry about?

3. Are there sample code using Cholesky and/or analytical form for doing fitting and predicting together over multivariate inputs?

Thanks,

Andy

[Bob Carpenter]

> On Apr 27, 2015, at 7:31 AM, andy...@cubist-asia.com wrote:
>
> Hi Linas/Bob,
>
> This thread has been extremely helpful for me to get up to speed on fitting and predicting GP "in one go". Thanks very much guys!
>
> Based on, FitPredict, I have tried to add multivariate input support and automatic relevance determination (Ch 14.2, 14.3, doc: 2.6.2)
>
> 1. As per suggestion in the manual, I have made the following changes:
>
> parameters {
> ....
> vector<lower=0>[D] rho;
> .....
> }
>
> model {
> ....
> matrix [N, N] Sigma;
> for (i in 1 : (N - 1)) {
> for (j in i : N) {
> vector[D] tmp_exp_term;
> tmp_exp_term <- rho .* (x[i] - x[j]);
>
> Sigma [i, j] <- eta_sq * exp (-dot_self(tmp_exp_term));
> // exp (-rho_sq * pow (x [i] - x [j], 2));
> Sigma [j, i] <- Sigma [i, j];
> }
> }
> ....
> ....
> rho ~ cauchy (0, 5);
> }
>
> I want to ask, from efficiency perspective, are there much better ways to do it without moving to Cholesky/analytical form?

Nope, not really.

> 2. Attached is a traceplot of selected paramters i have run, with N = 20, D=3. While the predicted y2 have converged quite well, trace of eta_sq and rho seems to have done less well. Is this something I should worry about?

Yes, there's a problem in that the rho and eta are going to be highly
correlated because of the basic form. So they're only weakly identified
by the priors, so you're going to need to run enough iterations to get
reasonable mixing. I don't know any reparameterizations that would help here.

> 3. Are there sample code using Cholesky and/or analytical form for doing fitting and predicting together over multivariate inputs?

You still need to build up the covariance matrix before
factoring it, at which point using the Cholesky factor should
be easy.

- Bob

[andy...@cubist-asia.com]

ok, thanks.

 

> 2. Attached is a traceplot of selected paramters i have run, with N = 20, D=3. While the predicted y2 have converged quite well, trace of eta_sq and rho seems to have done less well. Is this something I should worry about?


Yes, there's a problem in that the rho and eta are going to be highly
correlated because of the basic form.  So they're only weakly identified
by the priors, so you're going to need to run enough iterations to get
reasonable mixing.  I don't know any reparameterizations that would help here.

I tried to run it for different iter setting. But, it seems that sd of rho would only get worse if i increase the iter from 2k to 100k.

I am attaching the code, sample data I ran on, summary of the fit objects and traceplots. Again, the predicted y seems to be stable though.

One thing to note: the X1/predictors are highly correlated - output from R cor function:

          [,1]      [,2]      [,3]
[1,] 1.0000000 0.9813094 0.9459689
[2,] 0.9813094 1.0000000 0.9763934
[3,] 0.9459689 0.9763934 1.0000000

Would this be an issue? Or, I need a tighter prior on rho?


 

> 3. Are there sample code using Cholesky and/or analytical form for doing fitting and predicting together over multivariate inputs?

You still need to build up the covariance matrix before
factoring it, at which point using the Cholesky factor should
be easy.

Thanks. Will work through them.
 
Andy

- Bob

[Bob Carpenter]

> On May 6, 2015, at 3:53 AM, andy...@cubist-asia.com wrote:
>
>
>
> On Tuesday, May 5, 2015 at 4:43:44 AM UTC+8, Bob Carpenter wrote:

> ...

> I tried to run it for different iter setting. But, it seems that sd of rho would only get worse if i increase the iter from 2k to 100k.

I think the goal should be to get the right posterior sd,
so I'm not sure what you mean by worse --- is this versus
simulated data?

> I am attaching the code, sample data I ran on, summary of the fit objects and traceplots. Again, the predicted y seems to be stable though.

That happens when there is a range of parameter space that produces
similar outputs --- the predictions will have higher n_eff estimates
and better apparent convergence than their basis in coefficients, etc.

>
> One thing to note: the X1/predictors are highly correlated - output from R cor function:
>
> [,1] [,2] [,3]
> [1,] 1.0000000 0.9813094 0.9459689
> [2,] 0.9813094 1.0000000 0.9763934
> [3,] 0.9459689 0.9763934 1.0000000
>
> Would this be an issue? Or, I need a tighter prior on rho?

This can cause problems in fitting. The general thing to do is simulate
data like what you have and see if you can recover the correct values of
the parameters. Given that you're also trying to do prediction of new
y, that should be a reasonable way to evaluate whether what you're
trying to do will work.

- Bob

[andy...@cubist-asia.com]

On Wednesday, May 6, 2015 at 9:42:20 PM UTC+8, Bob Carpenter wrote:

> On May 6, 2015, at 3:53 AM, andy...@cubist-asia.com wrote:
>
>
>
> On Tuesday, May 5, 2015 at 4:43:44 AM UTC+8, Bob Carpenter wrote:
> ...
> I tried to run it for different iter setting. But, it seems that sd of rho would only get worse if i increase the iter from 2k to 100k.

I think the goal should be to get the right posterior sd,
so I'm not sure what you mean by worse --- is this versus
simulated data?

Actually, the predicted y_tilde matches quite well with the y.

But, my issue is more on rho - particularly seems the sd is quite unstable - at iter=2000, rho[3] has a mean and sd of 5.77 and 45.99 respectively. When i increase the iter to 100000, rho1, rho2, rho3 have (mena, sd) of (5.34, 140.49), (6.17, 140.89), and (6.65, 112.82) respectively.

Is this something I can fix through better model/prior specification or is it something I have to live with given the (highly) correlated predictors?

Many thanks,

Andy



 

[Bob Carpenter]

> On May 11, 2015, at 4:11 AM, andy...@cubist-asia.com wrote:
>
>
> On Wednesday, May 6, 2015 at 9:42:20 PM UTC+8, Bob Carpenter wrote:
>
> > On May 6, 2015, at 3:53 AM, andy...@cubist-asia.com wrote:
> >
> >
> >
> > On Tuesday, May 5, 2015 at 4:43:44 AM UTC+8, Bob Carpenter wrote:
> > ...
> > I tried to run it for different iter setting. But, it seems that sd of rho would only get worse if i increase the iter from 2k to 100k.
>
> I think the goal should be to get the right posterior sd,
> so I'm not sure what you mean by worse --- is this versus
> simulated data?
>
> Actually, the predicted y_tilde matches quite well with the y.
>
> But, my issue is more on rho - particularly seems the sd is quite unstable - at iter=2000, rho[3] has a mean and sd of 5.77 and 45.99 respectively. When i increase the iter to 100000, rho1, rho2, rho3 have (mena, sd) of (5.34, 140.49), (6.17, 140.89), and (6.65, 112.82) respectively.
>
> Is this something I can fix through better model/prior specification or is it something I have to live with given the (highly) correlated predictors?

You need to run until the estimates for the posterior means of rho1, rho2, and rho3
have converged. What's the r-hat and n_eff at 2000 (post-warmup!) iterations?

There are some densities, like the Cauchy, that don't have
finite mean and variance. It's unlikely to come up in a posterior
unless you have very fat-tailed priors and very little data.

- Bob

[Aki Vehtari]

There are some densities, like the Cauchy, that don't have
finite mean and variance.  It's unlikely to come up in a posterior
unless you have very fat-tailed priors and very little data.

Note that eta and rho are weakly identifiable as they are tied together. See 

 - Section 8.4 in Vanhatalo, J., Pietiläinen, V., & Vehtari, A. (2010). Approximate inference for disease mapping with sparse Gaussian processes. Statistics in medicine,29(15), 1580-1607.

and

 - Zhang H. Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. Journal of the American Statistical

Association 2004; 99(465):250--261. 

 - Diggle, P. J., Tawn, J. A., and Moyeed, R. A. (1998). Model-based geostatistics. Journal of the Royal Statistical Society. Series C (Applied Statistics), 47(3):299–350.

This may make the posterior difficult to sample unless good (weakly) informative priors are used.

Aki

[Andy Choi]

On Mon, May 11, 2015 at 11:49 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> On May 11, 2015, at 4:11 AM, andy...@cubist-asia.com wrote:
>
>
> On Wednesday, May 6, 2015 at 9:42:20 PM UTC+8, Bob Carpenter wrote:
>
> > On May 6, 2015, at 3:53 AM, andy...@cubist-asia.com wrote:
> >
> >
> >
> > On Tuesday, May 5, 2015 at 4:43:44 AM UTC+8, Bob Carpenter wrote:
> > ...
> > I tried to run it for different iter setting. But, it seems that sd of rho would only get worse if i increase the iter from 2k to 100k.
>
> I think the goal should be to get the right posterior sd,
> so I'm not sure what you mean by worse --- is this versus
> simulated data?
>
> Actually, the predicted y_tilde matches quite well with the y.
>
> But, my issue is more on rho - particularly seems the sd is quite unstable - at iter=2000, rho[3] has a mean and sd of 5.77 and 45.99 respectively. When i increase the iter to 100000, rho1, rho2, rho3 have (mena, sd) of (5.34, 140.49), (6.17, 140.89), and (6.65, 112.82) respectively.
>
> Is this something I can fix through better model/prior specification or is it something I have to live with given the (highly) correlated predictors?


You need to run until the estimates for the posterior means of rho1, rho2, and rho3
have converged.  What's the r-hat and n_eff at 2000 (post-warmup!) iterations?

Across all the runs, RHat are 1.

for iterations(post-warm-up) for 2000:n_eff are 2238, 3297, and 5221 respectively
for iterations(post-warm-up) for 50000:n_eff are 167077, 35449, and 71540 respectively

 

There are some densities, like the Cauchy, that don't have
finite mean and variance.  It's unlikely to come up in a posterior
unless you have very fat-tailed priors and very little data.

- Bob

--
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/Xt7JoUuOUkc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Andy Choi]

Hi Aki,

Many thanks. What kind of prior would you suggest for rho (in a multivariate setup)?

I have changed the prior from cauchy (0, 5) to normal (0, 5).

The mixing seems to be much better. But, I do note that I can move the prior from (0, 5) to (0, 100), the mean and sd of rho also tends to up.

Thanks,

Andy

--

[Aki Vehtari]

You didn't mention anything about your data or what prior information you have.

I usually scale x to have mean 0 and std 1. Then as a weakly informative prior I would use

 l ~ t(4,0,1)

where l=1/rho is the lengthscale. This prior is quite flat for l<1 and for l>1 the behavior of the GP does not change much (it's already close to linear), so it's unlikely that larger values are needed, but having the thick tail still provides the possibility that the data overrules the prior.

For eta I would use prior which says that eta_sq is likely to be smaller than the total variance of the data.  I usually scale y to have mean 0 and std 1, and thus I can use prior favoring eta<1. Note that, the total variance of the predicted data is eta_sq+sigma_sq, so you could also define a prior on eta_sq+sigma_sq and prior on ratio eta_sq/sigma_sq which is the proportion of the explained variance.

Aki

[Andy Choi]

Hi Aki,

Many thanks for the note.

Regarding the data:

y: normally distributed with mean of 0 and typical sd in the range of 0.1 to 0.5

x:
- each of the x is "assumed" to be normally distributed with mean 0 and sd related to the sd of y

- for those x(s) that are correlated/with explanatory power on y, they would tend to be highly correlated

functional form:

- my prior expectation is that y should be a linear function with zero intercept of x for  "small values" of x (most of the observations) but is bounded by an unknown asymptote values on both positive and negative ends. The intuition is that y will react to x linearly, but only up to a certain extend. Beyond that, an incremental change would have increasingly small impact on y.

Data size:

- number of observations(N): typically, 10 - 15

- number of x(D): can be up to 15

Here is a sample data set that I'm working off with:

$D
[1] 9

$N
[1] 10

$y
 [1] -0.030992628  0.148615099  0.094613001  0.003396619  0.111006356
 [6]  0.064603953 -0.007574454 -0.040076375 -0.128729695 -0.140480064

$x
$x[[1]]
[1] 0.007355116 0.012240551 0.029899214 0.021653077 0.034853025 0.028182862
[7] 0.041517091 0.017900408 0.005338227

$x[[2]]
[1] 0.07218964 0.05708903 0.06075094 0.06596872 0.05186838 0.06689468 0.05458556
[8] 0.05867080 0.03884169

$x[[3]]
[1] 0.03756785 0.01882525 0.04015725 0.03613047 0.04031964 0.03848678 0.03560452
[8] 0.03910110 0.02000822

$x[[4]]
[1] 0.06377052 0.05141590 0.05266848 0.04119564 0.04427342 0.04704632 0.05772354
[8] 0.06540720 0.03819622

$x[[5]]
[1] 0.05266989 0.04917466 0.04651094 0.03475764 0.04102571 0.04871280 0.05310855
[8] 0.04475245 0.02996702

$x[[6]]
[1] -0.005838582 -0.011384698 -0.014630407  0.024398174  0.003027790
[6] -0.010331219 -0.012495347 -0.001278804 -0.020244412

$x[[7]]
[1] -0.02012468 -0.01328256 -0.02256433 -0.02073496 -0.01436617 -0.01971798
[7] -0.01512621 -0.01654797 -0.01699086

$x[[8]]
[1] -0.04089289 -0.02621516 -0.04307371 -0.03965620 -0.04275576 -0.04397031
[7] -0.04919355 -0.04317085 -0.02979750

$x[[9]]
[1] -0.06954436 -0.05982779 -0.08587940 -0.11138529 -0.09631904 -0.08866991
[7] -0.06123835 -0.07207325 -0.03507037

$x[[10]]
[1] -0.10493675 -0.08487418 -0.13804823 -0.12619537 -0.12097182 -0.12688081
[7] -0.09680075 -0.10444238 -0.04869718

Many thanks,

Andy

[Aki Vehtari]

functional form:

- my prior expectation is that y should be a linear function with zero intercept of x for  "small values" of x (most of the observations) but is bounded by an unknown asymptote values on both positive and negative ends. The intuition is that y will react to x linearly, but only up to a certain extend. Beyond that, an incremental change would have increasingly small impact on y.

I assume you mean "close to linear" and not "linear". If you assume saturation, you should use some sigmoid function possibly combined with a GP, instead of just GP with a prior mean zero and a stationary covariance function.

 

Data size:

- number of observations(N): typically, 10 - 15

- number of x(D): can be up to 15

The latest model I found in this thread was t_gp_fp.stan. It has a GP prior with squared exponential covariance function allowing nonlinearities and full interactions. If rho is large for all dimensions, then you have too much flexibility compared to your small N. Since your x's are highly correlated, it would be useful use factor analysis distance (projecting x's to a lower dimensional manifold). Alternatively you need to have a strongly sparsifying prior for rho and hope that only one (or two) dimensions will be used. 

 Aki

[Andy Choi]

Hi Aki,

Many thanks.

1. I was meant to say "close to linear".

2. Yes, that the GP model I'm working right now. In terms of the rho's, will explore the approaches you suggested. But, just thinking out loud, would it work if I assume some kind of multivariate prior over the rho's?

[Aki Vehtari]

I checked the data you posted. There's not many observations, and it seems to be quite noisy. You mentioned that you assume saturation. I think you should start by fitting sigmoidal function and maybe with Student's t distribution first. It's likely that GP does not provide any improvement over that...

Aki

[Andy Choi]

Hi Aki,

Many thanks for the note.

Yes, I did more experiments and came to the same conclusion - I didn't get much (if any) improvement with GP over student-t. So, will prob go back to that for now.