Scaled Beta Prior

[Colin]

Dear Stan Group,

I am having trouble identifying the range parameter in a GP model with a Matern covariance function.  Priors such as a Cauchy or other diffuse priors on phi result in draws that are way outside the spatial range of the data.  I was thinking a scaled beta prior on phi (the range) parameter, could be a good option.  I was wondering if there was an option to do a scaled beta prior in Rstan or if there were other suggestions for alternative priors that would respect the spatial range of the data.

Thank you

[Aki Vehtari]

Note that we have moved the user forum to http://discourse.mc-stan.org/

We are just updating the GP section in the manual and adding prior recommendations https://github.com/stan-dev/stan/pull/2301. In the case you mention, our (default) recommendation is generalized inverse Gaussian distribution. There's not yet a builtin density function for that, although it will be added soon. 

Aki

[Colin]

Thanks, Aki.  I will post there moving forward.   I will look into the generalized inverse Gaussian distribution, too.  

In terms of the scaled beta distribution, would that be something I should just define as a new distribution in Stan and then add it as a prior?  I know Stan advocates using priors without bounded support, but in the case of a spatial model, I think there is a justification for using a prior that does not put mass on a range parameter that is beyond the maximum distance in the data.  Would you agree?

Best,

Colin

[Aki Vehtari]

Here's the Stan function for the generalized inverse Gaussian

functions {
  real generalized_inverse_gaussian_lpdf(real x, int p, real a, real b) {
      real lp =
      p * 0.5 * log(a / b)
      - log(2 * modified_bessel_second_kind(p, sqrt(a * b)))
      + (p - 1) * log(x) - (a * x + b / x) * 0.5;
      return lp;
  }
}

You can use it set a prior in the usual way

x ~ generalized_inverse_gaussian(p, a, b);

If the range parameter is the same as the maximum distance in the data then the correlation between the most distant points is still less than 1 and thus it is reasonable to "put mass on a range parameter that is beyond the maximum distance in the data". I'm not aware of cases where we would have just prior information that strict upper bound would be justified. The problem with a strict upper bound is that if you guessed wrong then no amount of data can overcome that. Bounds are ok, when the values outside the bounds are invalid (e.g. lower=0 for the range parameter, or lower=0,upper=1 for probabilities).

Aki

[Colin Lewis-Beck]

Thanks again, Aki,

I do see your point about the range parameter.   I will give this prior a shot.  The data set I have is pretty small (about 200 obs), so the posteriors for phi were extremely sensitive to the priors, and if the prior on phi was too diffuse the posterior was heavily right skewed with extremely large values.  

Below is the model I am fitting with a Matern covariance structure.   I am fixing the smoothness parameter to help simplify the model.  One idea I have to help the MCMC converge is to set the mean parameter of the phi prior equal to a value approximately equal to the range of the empirical variogram--sort of an empirical Bayes approach.   Do you think this would be a reasonable choice?

model {

   matrix[L,L] H; // Same Spatial Paramters for Each Year

   

   for (i in 1:(L-1)) {

     for (j in (i+1):L) {

       H[i,j] = tau2 * exp(-d[i,j]/ phi);

       H[j,i] = H[i,j]; // do not calculate twice, utilize symmetry

     }

   }

   for (k in 1:L) {

     H[k, k] = sigma2 + tau2; //diagonal of covariance matrix

   }

  

   for (i in 1:T){

   Y[,i] ~ multi_normal(alpha[i] + X[,i]*beta[i], H);

}

   

   //priors 

   //inproper priors on coefficients

   sigma ~ cauchy(0,1);

   tau ~ cauchy(0,1);

   phi ~ cauchy(0,5);

}

--
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/21qPDUWd9Fg/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.

[Aki Vehtari]

Please post all new questions to http://discourse.mc-stan.org/

Now you get only my answers, while in Discourse there would be many more reading these.

The next manual version will discuss sensitivity of covariance function priors in more detail. 

Cauchy is especially bad choice for phi.

Are you sure you want to use Matern with nu=1/2? That's far for from smooth. With 200 obs, nu will be weakly identifiable. I recommend testing with nu=3/2 or nu=5/2, too,

Setting the prior using empirical variogram is double use of data, and should be used only if the prior is wide enough so that it's still weakly informative.

You could also improve the convergence by using the non-centered parameterization.

Aki

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

[Colin Lewis-Beck]

Thanks, Aki.  I will post all future questions to the new board.   One final question if I may: when is the new manual planned for release?  I'd be interested in seeing what the default recommendations are for non-informative priors in the case of the GIG distribution.

Best,

Colin

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

[Bob Carpenter]

You don't need to assign---you can just return directly, so that'd
be

functions {
real generalized_inverse_gaussian_lpdf(real x, int p, real a, real b) {

return p * 0.5 * log(a / b)

- log(2 * modified_bessel_second_kind(p, sqrt(a * b)))
+ (p - 1) * log(x)
- (a * x + b / x) * 0.5;
}
}

I like to line up the additive terms on the log scale like this
to see the structure of the kernel; I'd probably also change x to y
and move the 0.5 out front, but that's even more cosmetic.

This could be rewritten to be a tiny bit more efficient by consolidating
terms inside the multiply by 0.5, but that hides the structure.

If you knew constant terms didn't matter, you could drop

* the 2* of the Bessel function
* if all of {p, a, b} are constants, the
first and econd terms can be dropped,
* if {p, x} are constants, then the third term can be dropped, and
* if {x, a, b} are constants, the fourth term can be dropped.

Our internal functions all do this, while at the same time being
careful to return the full result when used explicitly as functions
so as not to mess up things like mixtures.

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

[Aki Vehtari]

Colin, we prefer weakly informative priors over non-informative priors (especially as truly non-informative priors do not exist).

If you meant some default values for GIG parameters, the new manual doesn't have anything on those yet (only the recommended shape, ie GIG)

Bob, I copied the code from the manual and the line up was made so that it doesn't overflow to the margins and out of the page.

We have a problem there that we can't follow the code style recommendations as the stancode environment is using quite large font size compared to the page width.

And you are right that assignment could be removed.

Aki

[Bob Carpenter]

Thanks for looking into the formatting.

The real problem is the long function names. There's
a tension between not wanting to abbreviate and wanting
short names that's amplified with these compound names
and required prefixes.

Maybe that function should've been called mod_bessel_2nd or
even mod_bessel2 and maybe the defined function should be
gen_inv_gauss_lpdf or gen_inv_gaussian_lpdf or
general_inv_gaussian_lpdf.

Also, the Lucida console font we use for the actual manual
is narrower than Courier. It's one of the reasons it's a good
terminal font.

- Bob