I’m trying to create a simple model assuming data are distributed according to a lognormal distribution, but I want to put a prior on the mean of the lognormal rather than directly on the parameters of the lognormal. I created the mean as a transformed parameter and when using the ~ notation for the prior on the mean I got the warning:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
My question is, does simply changing from the ~ notation to target += fix the problem? The warning goes away, but I’m not sure if what I’m doing in appropriate.
data{
int<lower=0> n;
real<lower=0> y[n];
}
parameters{
real mu;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0> mean;
mean = exp(mu + sigma^2/2);
}
model {
target += normal_lpdf(mean | 10, 2); // as opposed to mean ~ normal(10,2);
y ~ lognormal(mu, sigma);
}
I’m not sure what’s Stan’s behavior when renaming existing functions, but since mean is a function that computes the mean of an array and is on the left hand side of the ~, do you still get the error if you name the variable something else?
Switching the notation doesn’t change anything. In both cases, you are incrementing the target by either normal_lpdf(mean | 10, 2) or normal_lupdf(mean | 10, 2). If your intention is that the prior pushforward distribution for mean should be a normal distribution, then you have two options. Either you can follow your current notation and include a Jacobian adjustment (see Jacobian adjustments explained simply), or you can declare mean as a parameter and derive mu as the transformed parameter, in which case no jacobian adjustment will be necessary.
Thank you, that’s exactly what I was looking for. I wasn’t familiar with this, and in looking this stuff up I felt like I was only seeing these transformations in the context of the response, so wasn’t sure if I needed to do anything given I was transforming a parameter and not my response.
I think your article was very helpful, but just to ensure I got what you were stating I have these two implementations to ensure that the prior pushforward distribution for mean is a normal distribution.
Option 1: Include Jacobian adjustment (adjusting with respect to mu and sigma since my transformation is nonlinear in both).
data{
int<lower=0> n;
real<lower=0> y[n];
}
parameters{
real mu;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0> mean;
mean = exp(mu + sigma^2/2);
}
model {
mean ~ normal_lpdf(10, 2);
target += log( exp(mu + sigma^2/2) ); //adjustment, partial with respect to mu
target += log( exp(mu + sigma^2/2) * sigma ); //adjustment, partial with respect to sigma
y ~ lognormal(mu, sigma);
}
Option 2: Declare mean as a parameter and derive mu as the transformed parameter
data{
int<lower=0> n;
real<lower=0> y[n];
}
parameters{
real mean;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0> mu;
mu = log(mean) - sigma^2/2;
}
model {
mean ~ normal_lpdf(10, 2);
y ~ lognormal(mu, sigma);
}
I simulated data from these and the first implementation looks correct, the prior on mean appears to be distributed as I wanted to imply.
The second option doesn’t seem to be correct. I assume I would need to do a Jacobian adjustment there but can’t seem to get that correct.
I pretty sure trying to find a Jacobian is fruitless as the transformation from mu and sigma to the mean is not a 1-to-1 function.
The approach I use in my paper is to define a loss function to induce a posterior which has the same density as your “expert opinion”. In this case it is best to treat E_mean (the mean of the log normal distribution) as a parameter with mu a function of E_mean and sigma (which should be uniform). You can either put a improper uniform prior on E_mean and update with a “loss function” or set E_mean to have the density of your expert opinion (a Normal distribution), in which case there is no-loss function required.
Clearly the second approach is more intuitive, however, the loss framework is more general as it can include information when the observable outcome is not a straightforward function of the other parameters.
Loss approach
model_test2<- "
parameters{
real<lower=0> E_mean; //don't use the name "mean" as it is a function
real<lower=0> sigma;
}
transformed parameters{
real mu;
mu = log(E_mean) - sigma^2/2;
}
model {
// sigma can even be an improper uniform (should be non-informative)
sigma ~ uniform(0,10);
target += normal_lpdf(E_mean | 10, 2); //(negative) Loss function
}
"
Hierarchical approach
model_test<- "
parameters{
real<lower=0> E_mean;
real<lower=0> sigma;
}
transformed parameters{
real mu;
mu = log(E_mean) - sigma^2/2;
}
model {
sigma ~uniform(0,10);
E_mean ~ normal(10, 2);
}
"
Thank you, I didn’t think about the use of mean in writing this post. I didn’t actually use mean in production so I didn’t get any error or anything related to that when trying to compile and execute my actual code. But that’s a good observation I didn’t consider in preparing my post.
MCMT_.5_1.pdf (865.1 KB)
Thanks. This was helpful and appears to provide the solutions I was looking for. I uploaded a document showing 5 different implementations I tried and the two you suggested appear to replicate the prior distribution I desired whereas the other implementations do not appear to. In the “prior modeling” of the loss and hierarchical approaches you suggested I get lots of divergent transitions but I don’t think I care about that. I really just wanted to make sure the model I specified matched what I desired and thought I was specifying, and sampling the priors suggest that I do get what I desired. When I incorporated some simulated data to include the likelihood contribution to the posterior I don’t get the divergent transitions.