Need some help handling a transformation (inverses?) in a hierarchical bayesian model

[David Siegel]

Hey folks

I'm new to bayesian modeling with Stan, and after reading the forum and relevant parts of the manual (I hope) I'm not totally sure if I'm properly implementing a transformation of a variable.  I've copy-pasted my full code below.

The fundamental issue is that I have a prior on a transformed variable, but the inverse transformation cannot be written down with pen and paper (although it's probably invertible).

To give a simple illustration of the problem, imagine I have something like:

Simple Pseudocode:

data {

    int<lower=0> N;

    int<lower=1> K;

    vector[N] y; 

    matrix[N,K] x;

}

parameters {

    vector[K] beta;

    real a0;

}

model {

    f(a0) ~ uniform(log(1),log(N)); // f is an invertible transform from reals to reals, but it's inverse can't be written down

    target += log(f'(a0)); // The jacobian adjustment

    target += g(a0,betas,x,y);  // g is some complicated expression -- the point is that I need to use a0 here

}

Is it ok to continue to use a0 in the third line of the model, or do I need to write things in terms of the inverse of f?  If it's the latter, then the problem is that it's impossible for me to write down the inverse of f with pen and paper -- is there a workaround?

My code below is a little more complicated than this, but I think this captures the basic idea.  In the above example, "f" would be "log_sum_pis" below.

Thank you for your help!

David

Actual Code:

model_code = """    

data {

    int<lower=0> N;

    int<lower=1> K;

    vector[N] y;

    matrix[N,K] x;

    real<lower=0> lambda;

    vector[K] exp_term;

}

parameters {

    vector[K] beta;

    real a0;

}

transformed parameters {

    real<lower=0> squared_error;

    real<lower=0> lambda_delta;

    vector[K] pis;

    vector[K] pi2exp;

    real log_sum_pis;

    squared_error = dot_self(y - x * beta) / N;

    lambda_delta = lambda * 100;

    for (i in 1:K)

        pis[i] = 1/(1+exp(-a0-exp_term[i]));

    for (i in 1:K)

        pi2exp[i] = exp(-exp_term[i])/(1+exp(-a0-exp_term[i]))/(1+exp(-a0-exp_term[i]));

    log_sum_pis = log(sum(pis));

}

model {

    log_sum_pis ~ uniform(log(1),log(N));

    target += log(sum(pi2exp)/sum(pis)); //The Jacobian adjustment for the above line

    target += -squared_error;

    for (k in 1:K)

        target += log_sum_exp(log(pis[k]*lambda/2) - lambda * fabs(beta[k]),log((1-pis[k])*lambda_delta/2) - lambda_delta * fabs(beta[k])); //Note that this line uses "pis" and not "log_sum_pis".

}

"""

[Andrew Gelman]

HI, I think you want to avoid this uniform distribution.

A

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

[Bob Carpenter]

Yes, that's fine technically.

Just because f(a0) is uniform on some interval, doesn't mean
a0, the parameter of interest is. I'm guessing f() is just
a convenient transform to define the prior density you want.

If a0 winds up being constrained, and the constraint isn't
a hard one such as covariance matrices being positive definite
or probabilities being in [0, 1], then that can lead to computational
and statistical problems at the boundaries.

- Bob

> On Feb 27, 2017, at 11:04 PM, David Siegel <davida...@gmail.com> wrote:
>

[David Siegel]

Ah thank you.

Yes, I think you've got it.  I'm really trying to define the prior distribution of a0, but f() is the thing I actually know the distribution of.

Thanks for the heads up.  And thank you for the response!

David