Lognormal vs taking the exponential of the parameter

Modeling
1

Given a parameter that I expect to follow a lognormal distribution (e.g., because I know the value must be greater than 0), I find it more intuitive to use a transformation. So, instead of doing the following where I need to use a complicated formula to convert the values of mu and sigma so that I can obtain the mean value of the distribution:

parameters {
    real mu;
    real sigma;
    real gain;
}

model {
    gain ~ lognormal(mu, sigma);
}

Can I just do the following?

parameters {
    real mu;
    real sigma;
    real t_gain;
}

transformed parameters {
    real gain;
    gain = exp(t_gain);
}

model {
    t_gain ~ normal(mu, sigma);
}

This is just a snippet of a much larger model, but I’m wondering if my approach is reasonable or if there’s a reason to favor the added conceptual complexity of the lognormal approach? Is the lognormal more likely to converge or sample better?

1 Like
2

I think you are mixing up.

Log normal is

log(y) ~ normal(…, …)

If this is the case, you have to apply the Jacobian. Long story short, use

y ~ lognormal(…, …);

applying the right jacobian, no difference. Just use the lognormal ditribution. They have worked it out for you.

3

I’m a little confused. If we are to do:

exp(log(y)) ~ exp(normal(..., ...))

This gives us:

y ~ exp(normal(..., ...))

Which is what I’ve proposed in my original post (except that the transformation is performed in the transformed parameters block).

Are they mathematically equivalent? Do I need to perform a Jacobian correction regardless of whether I do log(y) ~ normal(..., ...) or y ~ exp(normal(..., ...))?

4

Those two models are equivalent. Though in the first model, you would want to write
real <lower=0> gain; rather than real gain;

The second model will likely be more numerically stable, since under the hood when you do the greater than zero constraint, it does the exponential transformation, and then takes the log of x in the lognormal distribution statement.

5

To confirm, you feel that the second model (i.e., the one I pasted below) is the better one to use? No Jacobian correction is needed, right?

parameters {
    real mu;
    real sigma;
    real t_gain;
}

transformed parameters {
    real gain;
    gain = exp(t_gain);
}

model {
    t_gain ~ normal(mu, sigma);
}
6

Correct. The Jacobian corrected model would be:

parameters {
    real mu;
    real sigma;
    real<lower=0> gain;
}

model {
    target += -log(gain);
    log(gain) ~ normal(mu, sigma);
}

Which you can see is identical to the model with gain ~ lognormal(mu,sigma).

1 Like
7

The function “normal” does produce log probability values (log dnorm in R) not actual normally distributed values (rnorm in R ).

Have a read to the function manuAl of stan. You will quickly as we to your ouw doubts.

8

@stemangiola: Are you fine with the approach proposed by @aaronjg? He seems to think it’s a bit more computationally stable. I’m running it now, but my model is very slow so I won’t know for some time whether it works.

9

Yes

10

Hi, a follow on question:

I have a lognormal model that requires addition of two components. Is there any adjustment needed for the following?

data {
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
}
parameters {
  real log_beta;
  real<lower=0> sigma_log_e;
}
model {
  vector[N] mu = x + exp(log_beta);
  y ~ lognormal(log(mu), sigma_log_e);
}

*Disclaimer: not a real model.

11

You mean this to be a linear model?

y ~ lognormal(x * beta, sigma);

Otherwise doesn’t make much sense (unless I am missing something)

12

Sorry, I was trying to be brief. More like this:

parameters {
  real log_beta;
  real<lower=0> sigma_log_e;
  vector<lower=0, upper=1> [N] prop;
}
model {
  vector[N] mu = prop .* x + (1 - prop) * exp(log_beta);
  y ~ lognormal(log(mu), sigma_log_e);
}

Again, not a real model but an attempt to communicate that addition is necessary.

I’m not across fully my Jacobian adjustments, but does it count as a change of variables if I exponentiate and transform back?

13

Following the simple roule that any transformation to the right side of the ~ does not requite Jacobian adjustments, your case does not require it.

1 Like
14

No adjustment needed. Though it may be more numerically stable to transform x to log_x in the transformed data block and use log_mix(prop,log_x[i],log_beta[i]). On the other hand, it won’t be vectorized so may be slower. If I were you, I’d keep it as is for now, but if you start running into divergent transitions, you should consider using log_mix. You might even want to roll your own version of it, and handle the translation of prop_raw on the unconstrained space to prop on the constrained space and then calculate the appropriate values for log(prop) and log(-prop) using prop_raw.

Anyway, probably not worth worrying about unless you see stability issues, as indicated by divergences.

1 Like