Loo: Adjust log-lik for non-linear transformation

General
1

I’m trying to compare two models in Stan using the loo package in R. Fitting the models is no problem. The only issue is making sure that I calculate the log-likelihood correctly in one of them. I would imagine that getting the log likelihood wrong would be an issue when comparing the models.

The baseline model is just a simple univariate normal distribution (call the original variable Y). The more sophisticated version to splits the model up into two parts, let’s call them X_1 and X_2. Both X_1 and X_2 are fit with univariate normal distributions, but the second part requires a non-linear transformation to get back to Y (e.g. Y=X_1-log(1+X_2/K), where the parentheses and K are above zero, also expressed in terms of X_2 we would have X_2=(exp(X_1-Y)-1)*K)

My question is: 1) do I need to make an adjustment to the log-likelihood to account for this non-linear transformation, and 2) if so, what type of adjustment would I need to use?

A more fleshed out example is below (though with a different naming convention, the temp variable in generated quantities is the , chg_ln_P lines up with Y, chg_ln_D lines up with X_1, and chg_r_e lines up with X_2).

Thanks.

data {
    int<lower=1> T;
    vector[T] chg_ln_P;
    vector[T] chg_ln_D;
    vector[T + 1] r_e;
    real<lower=0> r_f;
    real<lower=0> r_e_star;
    real<lower=0> g;
}
transformed data {
    vector[T] chg_r_e;
    for (i in 1:T) {
        chg_r_e[i] = r_e[i + 1] - r_e[i];
    }
}
parameters {
    real mu_chg_ln_D;
    real alpha_chg_r_e;
    real beta_chg_r_e;
    real<lower=0> sd_chg_ln_D;
    real<lower=0> sd_chg_r_e;
}
model {
    mu_chg_ln_D ~ normal(0, 0.1 /12);
    alpha_chg_r_e ~ normal(0, 0.1 /12);
    beta_chg_r_e ~ normal(0, 1);
    sd_chg_ln_D ~ cauchy(0, 0.01);
    sd_chg_r_e ~ cauchy(0, 0.01);

    chg_ln_D ~ normal(mu_chg_ln_D, sd_chg_ln_D);
    chg_r_e ~ normal(alpha_chg_r_e + beta_chg_r_e * (r_e[1:T] - r_e_star), sd_chg_r_e);
}
generated quantities {
    vector[T] log_lik;
    for (i in 1:T) {
        log_lik[i] = normal_lpdf(chg_ln_D[i] | mu_chg_ln_D, sd_chg_ln_D);
        log_lik[i] = log_lik[i] + normal_lpdf(chg_r_e[i] | alpha_chg_r_e + beta_chg_r_e * (r_e[i] - r_e_star), sd_chg_r_e);
    }
}
2

I figured out the transformation that makes it work. The adjusted generated quantities is below.

Basically I needed to figure out the mu and sigma that work with exp(X_1-Y)~normal(mu, sigma) and then add in the transformation for exp(X_1-Y). So instead of a log-normal, this is kind of an exponential-normal. I don’t recall seeing that before, but it works for my case.

generated quantities {
    vector[T] log_lik;
    {
        vector[T] mu_exp;
        vector[T] sd_exp;
        mu_exp = 1 + (alpha_chg_r_e + beta_chg_r_e * (r_e[1:T] - r_e_star)) ./ (r_f + r_e[1:T] - g);
        sd_exp = sd_chg_r_e ./ (r_f + r_e[1:T] - g);
        for (i in 1:T) {
            log_lik[i] = normal_lpdf(chg_ln_D[i] | mu_chg_ln_D, sd_chg_ln_D);
            log_lik[i] = log_lik[i] + normal_lpdf(exp(chg_ln_D[i] - chg_ln_P[i]) | mu_exp[i], sd_exp[i]);
            //new?
            log_lik[i] = log_lik[i] + log(exp(chg_ln_D[i] - chg_ln_P[i]));
        }
    }
}
1 Like
3

There is an example when Jacobian needs to be included in Regression and Other Stories: Mesquite, and the book chapter has a bit more info.

1 Like
4

Cool. Thanks.