Using Gaussian Copulas for hierarchical priors

[Ben Lambert]

Hi,

This probably is a terrible idea, but I was wondering, is there any precedent for using Gaussian/other-Copulas to specify hierarchical priors?

The particular example I am considering is probably a model I wouldn't necessarily code in Stan, since the parameters are discrete. (Obviously this could be done by marginalising parameters, but it might be quite messy). However, I wanted to see whether the community 'Bayes-brain' had seen something like this before.

I have two parameters X1 and X2, which are both discrete, and I want to assign very different negative binomial marginal priors to each. However, I would like there to be a rather complex dependence between these parameters. One way of specifying this dependence would be to use a Copula. This would allow me the flexibility to set the marginal priors independently, whilst allowing for an arbitrary dependence between the priors.

I know this can be done nicely for likelihoods, but I was wondering whether anyone had seen anything like this for priors? I suspect that it may suffice to use simpler options for examples where there is a lot of data. However in a few cases I am considering there are relatively few data points, and I would want to incorporate this complex dependency structure that is quite well known a priori.

Best,

Ben

[Ben Goodrich]

On Friday, February 19, 2016 at 9:52:30 PM UTC-5, Ben Lambert wrote:

This probably is a terrible idea, but I was wondering, is there any precedent for using Gaussian/other-Copulas to specify hierarchical priors?

 
I don't think it is a terrible idea, but a Gaussian copula doesn't permit a particularly complex dependency structure.

Ben

[andre.p...@googlemail.com]

Ben,

take a look at page 5 in:

http://support.sas.com/resources/papers/proceedings11/355-2011.pdf

I'm not exactly into your problem. The paper is useful in about how to

calculate the copulas for discrete distributions. It works very well.

Although calculation in STAN due to the nature of the copulas is not the 

fastest. 

Andre

[Ben Lambert]

Thanks both for your answers, and Andre for the reference. Still not sure on whether this is the best solution, but most multidimensional discrete distributions that I can find either aren't part of any MCMC packages, or just aren't sufficiently flexible. Any ideas still welcome! Best, Ben

[andre.p...@googlemail.com]

Ben,

the "best" solution is problem dependent. You can code the multidimensional distributions in STAN and compare their log-likelihoods, 

and/or additionally the number of parameters they use (->WAIC).

Andre

[John Hall]

What if you're not particularly interested in a complex dependency structure? Instead, you're just interested in the correlation between parameters with different uni-variate distributions. For instance,consider the case where one set of parameters are normally distributed and another set are log normally distributed and that they have some kind of correlation. One model I'm currently looking has something like this going on.

[Bob Carpenter]

Usually you work on the unconstrained scale and then map
back to the constrained scale.

So if you have y = (y[1], y[2]), where y[1] is unconstrained
and y[2] is constrained to be positive, then you model

(y_raw[1], y_raw[2])

as having a multivariate normal structure and define

y[1] = y_raw[1]

y[2] = exp(y_raw[2])

Same deal if some are probabilities (use inverse logit
instead of exp).

- 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[John Hall]

Yes, you're right. I'm pretty sure this is actually equivalent to using a Gaussian copula. It turns out that my data actually could probably be described with Gumbel copulas or something, but it's probably a small enough issue that I can ignore.

> 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/uGZAOemKt1s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[John Hall]

I know that the posterior distribution of the second set of parameters is log normal, but I might want to give them a half-normal prior (for the same reasons that it can be tricky to give a hierarchical variance a log normal prior). I'm thinking of using something like what's below to constrain y[2] to be bounded by 0. Would this have the right effect?

y ~ multi_normal(mu, sigma);

if (y[2] < 0)

  target += negative_infinity();

else

  target += -normal_lccdf(0 | mu[2], sqrt(sigma[2, 2]));

[Bob Carpenter]

Yes and no. It'd reject values below 0, but it's a
terrible way to code this.

Instead, declare y[2] with a lower bound of zero. If
the other y don't have constraints, paste together
a vector of unconstrained and constrained variables.

Anything else will kill sampling performance. Stan programs
should have finite density for any values of the parameters
that satisfy the declared constraints. Otherwise, they may
just devolve to simple random walks.

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

[John Hall]

You're talking about something like below.

parameters {
    vector[N] alpha;
    vector<lower=0>[N] beta;
    ...
}
transformed parameters {
    vector[N] gamma[2];
    ...
   
    gamma[1] = alpha;
    gamma[2] = beta;
    ...
}
model {
    gamma ~ multi_normal(mu_gamma, sigma_gamma);
    ...
}

The only problem I can see is that it becomes complicated to implement something like this with a multivariate reparameterization. Would something like the below work? In the more general case, I would be fitting a prior to mu and L_sigma, so delta[2]'s lower bounds would depend on those. That gets complicated.

parameters {
    vector[N] delta[2];
    ...
}
transformed parameters {
    vector[N] alpha;
    vector<lower=0>[N] beta;
    vector[N] gamma[2];
    ...
   
    gamma = mu_gamma + L_sigma_gamma * delta;
    alpha = gamma[1];
    beta = gamma[2];
    ...
}
model {
    delta ~ normal(0, 1);
    ...
}

> > 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/uGZAOemKt1s/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/uGZAOemKt1s/unsubscribe.

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

[Bob Carpenter]

> On Mar 20, 2017, at 11:32 AM, John Hall <john.mic...@gmail.com> wrote:
>
> You're talking about something like below.
>
> parameters {
> vector[N] alpha;
> vector<lower=0>[N] beta;
> ...
> }
> transformed parameters {
> vector[N] gamma[2];
> ...
>
> gamma[1] = alpha;
> gamma[2] = beta;

If you want the positive constrained thing to have
a truncated multivariate normal, then yes, but you
probably have the multi-normal backward, because presumably
it's the (alpha, beta) pairs that get a multinormal.

But the more common way to do this is to have the multi-normal
work on the unconstrained forms and then transform afterward
with exp() or inv_logit() to get positive or (0, 1) constrianed
variables:

vector[2] gamma_unc[N];

gamma_unc ~ multi_normal(...);

vector[2] gamma[N];
gamma[ , 1] = gamma_unc[ , 1];
gamma[ , 2] = exp(gamma_unc[ , 1]);

Then the multivariate reparameterization works as normal with
gamma_unc.

- Bob