Two simple equivalent (?) log-normal regrssion models -> different results!

[Jan]

Hello,

I have two representation of the same regression

y = a * x + mu

y - a * x ~ lognormal(mu, sigma);

1) works but gives Warning

2) is stable but give a slightly wrong result

1)------------------------------------------------

data{

  int N;

  vector[N] y;

  real mu;

  real sigma;

  vector<lower=0>[N] x;

}

parameters{

  real<lower=0, upper=1> a;

}

transformed parameters{

  vector[N] y_regressed_log;

  for(n in 1:N) y_regressed_log[n] = log(y[n]) + log1m((a*x[n]) ./ y[n]);

}

model{

y_regressed_log ~ normal(mu, sigma);

for(n in 1:N) target += -log(y[n]) - log1m((a*x[n]) ./ y[n]);

}

The following numerical problems occured the indicated number of times after warmup on chain 4
                                                                                            count
Exception thrown at line 19: normal_log: Random variable[2] is nan, but must not be nan!       28
Exception thrown at line 19: normal_log: Random variable[1998] is nan, but must not be nan!     5
Exception thrown at line 19: normal_log: Random variable[21] is nan, but must not be nan!       4

2)-----------------------------------------------------

data{

  int N;

  vector[N] y;

  real log_mu;

  real log_sigma;

  vector<lower=0>[N] log_x;

}

parameters{

  real log_a;

}

transformed parameters{

  vector[N] log_model;

  for(n in 1:N){

    vector[2] temp;

    temp[1] = log_a + log_x[n];

    temp[2] = log_mu;

    log_model[n] = log_sum_exp(temp);

  }

}

model{

y ~ lognormal(log_model, log_sigma);

}

-----------------------------------------------------------------------

Any ides?

Thanks a lot!

[Bob Carpenter]

> On Aug 5, 2016, at 8:30 AM, Jan <mangiol...@gmail.com> wrote:
>
> Hello,
>
> I have two representation of the same regression
>
> y = a * x + mu
>
> y - a * x ~ lognormal(mu, sigma);
>
> 1) works but gives Warning
> 2) is stable but give a slightly wrong result

I didn't follow along in detail, but I don't think your (2) is
a version of your (1), which would explain why they get different
answers. You need the Jacobian term to transform from normal
to lognormal parameterizations.

You can try to increase adapt_delta to avoid overflow or
underflow (which is where those warnings are coming from).
But they're no big deal.

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

[Jan]

On Sunday, August 7, 2016 at 10:57:43 AM UTC+10, Bob Carpenter wrote:

> On Aug 5, 2016, at 8:30 AM, Jan <mangiol...@gmail.com> wrote:
>
> Hello,
>
> I have two representation of the same regression
>
> y = a * x + mu
>
> y - a * x ~ lognormal(mu, sigma);
>
> 1) works but gives Warning
> 2) is stable but give a slightly wrong result

I didn't follow along in detail, but I don't think your (2) is
a version of your (1), which would explain why they get different
answers.  You need the Jacobian term to transform from normal
to lognormal parameterizations.

That's right

indeed if I log() a linear regression having lognormal shape the result is not even lognormal anymore. So the question becomes: what is the Jacobian of 

log(y - a * x) ~ normal(mu, sigma) 

 

is

-log(y[n]) - log1m((a*x[n]) ./ y[n])  ?

Thanks!

[Bob Carpenter]

Before trying to figure out a Jacobian, what are
you trying to do? It's odd to see a regression
on the left-hand-side of ~ because it won't arise
from a generative model. Shouldn't just the y be
on the left side and the a * x rolled into the mean
of the lognormal?

- Bob

[Jan]

I have a quite peculiar linear regression but with a lognormal/gamma noise model.

observed_ln = ln_1 * proportion_1 + ln_2 * (1 - proportion_1)

Please have a look to this tread https://groups.google.com/forum/#!category-topic/stan-users/general/jQKzA6vfgtU . Few questions are about solving the same problem in different ways

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

Hi Jan,

 

I have two representation of the same regression

y = a * x + mu

y - a * x ~ lognormal(mu, sigma);

1) works but gives Warning

2) is stable but give a slightly wrong result

Assuming y is distributed log-normal. (please specify more precisely),

then a*x + mu is on exp(a*x+mu) scale in 1), see:

https://en.wikipedia.org/wiki/Log-normal_distribution

and in 2) -a*x is on log-scale whereas only mu is on exp(scale).

Is it correct to say, you want to do a log-normal regression in STAN

with only using a normal distribution and asking for the jacobian transformation

to do so, because you need some parameter restrictions which cannot

be done by using STAN lognormal function?

Andre