Jacobian adjustment for parameters with parameter-based upper/lower bounds

Modeling
1

Hi all,

I’m trying to understand whether Jacobian adjustment is needed when a parameter is used as the upper and/or lower bound for other parameters.

For example,

parameters {
  vector<lower=0, upper=10>[k] a1;
  vector<lower=a1, upper=10>[k] a2;
  vector<lower=a2, upper=10>[k] a3;

  vector<lower=0, upper=10>[k] b1;
  vector<lower=b1, upper=10>[k] b2;
  vector<lower=b1 + b2, upper=10>[k] b3;
}

model {
  a1 ~ normal(0, 2);
  a2 ~ normal(0, 2);
  a3 ~ normal(0, 2);
  
  b1 ~ normal(0, 2);
  b2 ~ normal(0, 2);
  b3 ~ normal(0, 2);
}

I believe that when the upper and lower bounds are constants, as in a1 and b1, no additional Jacobian adjustment is needed.

My questions are:

  • Does Stan automatically handle the constraint transformation for a2 and a3, or do we need to add an explicit Jacobian adjustment in the model block?

  • Does the answer change when the upper and lower bounds are parameters, for example vector<lower=a1, upper=a3>[k] a2?

  • Does it make a difference if we include multiple parameters with parameter-dependent upper and lower bounds, as shown in the code for b3?

Any clarification would be greatly appreciated.

Thanks.

2

When using the built-in constraints in the parameters block the Jacobian adjustments are performed automatically.

For upper or lower bound constraints the Jacobian adjustment is non-linear and whether the bounds are parameters or not doesn’t matter. Typically the adjustment is desirable because then you can use the adjusted parameter on its constrained domain like any other parameter in your model. There are occasions when placing a prior on the unadjusted parameter(s) are preferred but these would be declared without bounds or other constraints. For example, take a standard deviation parameter sigma, I might want to place a skew_normal prior on log(sigma). I could declare sigma with a lower bound of 0 then take the log of it and place the prior on this but I would need a Jacobian adjustment. More efficiently I could declare an unbounded parameter log_sigma and place the prior directly on this. Subsequently, I would put exp(log_sigma) as the standard deviation into the normal density (let’s assume I’m modeling the s.d. of the normal) and I don’t need any Jacobian adjustment.

3

Thanks @spinkney

In my case, the parameters are regression coefficients. Let’s say a1, a2, and a3 are the raw, unconstrained parameters, and then in the transformed parameters block I construct b1, b2, and b3 using parameter-dependent upper and lower bounds.

For example, something like:

parameters {
  vector[k] a1;
  vector[k] a2;
  vector[k] a3;
}
transformed parameters {
  vector<lower=0, upper=a2>[k] b1;
  vector<lower=a1, upper=a3>[k] b2;
  vector<lower=a2, upper=10>[k] b3;
}
model {
  a1 ~ normal(0, 2);
  a2 ~ normal(0, 2);
  a3 ~ normal(0, 2);
}

In this setup, my priors are on the unconstrained a1, a2, and a3, not on the bounded b1, b2, and b3. In that case, would I still need to include an explicit Jacobian adjustment, if so, then could you please show how to add the Jacobian

My understanding is that if the priors are assigned to the unconstrained parameters, then the bounded parameters are just deterministic transformations, so the Jacobian would only be needed if I wanted to write the prior on the constrained scale instead. Is that correct?

I want b1, b2, and b3 to be ordered / constrained 0 < b1 < b2 < b3 < 10
Could you please confirm whether this transformed-parameters setup is correct in Stan?

Thanks

4

You can use the _jacobian functions to perform the adjustments as I do in Random walk with positive differences and missing data - #4 by spinkney. You can find more documentation at Variable Transformation Functions.

5

Thank you for the pointer — that’s very helpful. I’ll look into the these functions.

Regards

6

Yes.

Yes, it will depend on those parameters in this case.

It’s legal, but you should think through what the implied joint density is based on the Jacobian adjustments implied by the adjustments. Stan even lets you write things like this:

y ~ normal(0, 1);
y ~ normal(0, 1);

If you work through the algebra, the result is equivalent to

y ~ normal(0, 1 / sqrt(2));

What you posted is not a well-formed Stan model because b1, b2, and b3 are declared, but not defined. If you run it, the b values will all be not-a-number. So I’m not sure what the intent is here.

7

Thanks @Bob_Carpenter

I’m fitting a hierarchical piecewise spline model with subject-specific (random) change points. The outcome is a pain score measured repeatedly over time for each individual. The time variable is denoted C_1​, and the change points are the parameters b_t1​,...,b_t6,​, which are constrained to lie within certain bounds.

Below are my current parameters and transformed parameters blocks. The model code was generated using brms, which manually adds truncated-normal log densities to lprior for the bounded change-point parameters b_t1​,...,b_t6. Do I need to add any additional Jacobian adjustments for these bounded btb_tbt​ parameters? Key points about my setup:

  • The b_t1​,...,b_t6 parameters are declared directly as bounded parameters (with data- and parameter-dependent bounds), not as nonlinear transforms of other parameters.

  • brms is already assigning truncated-normal priors to them via explicit log-density terms in lprior.

  • I do not apply any further nonlinear transformation to b_t1​,...,b_t6 before using them in the likelihood.

Given this, is the current truncated-normal implementation sufficient, or is an extra Jacobian adjustment still required? I’d appreciate guidance on that as well.

Thanks

parameters {
  vector[K_cons] b_cons;  // regression coefficients
  vector[K_s1] b_s1;  // regression coefficients
  vector[K_s2] b_s2;  // regression coefficients
  vector[K_s3] b_s3;  // regression coefficients
  vector[K_s4] b_s4;  // regression coefficients
  vector[K_s5] b_s5;  // regression coefficients
  vector[K_s6] b_s6;  // regression coefficients
  vector[K_s7] b_s7;  // regression coefficients
  vector<lower=min(C_1),upper=max(C_1)>[K_t1] b_t1;  // regression coefficients
  vector<lower=min(C_1) + b_t1,upper=max(C_1)>[K_t2] b_t2;  // regression coefficients
  vector<lower=min(C_1) + b_t2,upper=max(C_1)>[K_t3] b_t3;  // regression coefficients
  vector<lower=min(C_1) + b_t3,upper=max(C_1)>[K_t4] b_t4;  // regression coefficients
  vector<lower=min(C_1) + b_t4,upper=max(C_1)>[K_t5] b_t5;  // regression coefficients
  vector<lower=min(C_1) + b_t5,upper=max(C_1)>[K_t6] b_t6;  // regression coefficients
  real<lower=0> sigma;  // dispersion parameter
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  matrix[M_1, N_1] z_1;  // standardized group-level effects
  cholesky_factor_corr[M_1] L_1;  // cholesky factor of correlation matrix
}
transformed parameters {
  matrix[N_1, M_1] r_1;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_1] r_1_cons_1;
  vector[N_1] r_1_s1_2;
  vector[N_1] r_1_s2_3;
  vector[N_1] r_1_s3_4;
  vector[N_1] r_1_s4_5;
  vector[N_1] r_1_s5_6;
  vector[N_1] r_1_s6_7;
  vector[N_1] r_1_s7_8;
  vector[N_1] r_1_t1_9;
  vector[N_1] r_1_t2_10;
  vector[N_1] r_1_t3_11;
  vector[N_1] r_1_t4_12;
  vector[N_1] r_1_t5_13;
  vector[N_1] r_1_t6_14;
  // prior contributions to the log posterior
  real lprior = 0;
  // compute actual group-level effects
  r_1 = scale_r_cor(z_1, sd_1, L_1);
  r_1_cons_1 = r_1[, 1];
  r_1_s1_2 = r_1[, 2];
  r_1_s2_3 = r_1[, 3];
  r_1_s3_4 = r_1[, 4];
  r_1_s4_5 = r_1[, 5];
  r_1_s5_6 = r_1[, 6];
  r_1_s6_7 = r_1[, 7];
  r_1_s7_8 = r_1[, 8];
  r_1_t1_9 = r_1[, 9];
  r_1_t2_10 = r_1[, 10];
  r_1_t3_11 = r_1[, 11];
  r_1_t4_12 = r_1[, 12];
  r_1_t5_13 = r_1[, 13];
  r_1_t6_14 = r_1[, 14];
  lprior += normal_lpdf(b_cons | 0, 1);
  lprior += normal_lpdf(b_s1 | 0, 5);
  lprior += normal_lpdf(b_s2 | 0, 5);
  lprior += normal_lpdf(b_s3 | 0, 5);
  lprior += normal_lpdf(b_s4 | 0, 5);
  lprior += normal_lpdf(b_s5 | 0, 5);
  lprior += normal_lpdf(b_s6 | 0, 5);
  lprior += normal_lpdf(b_s7 | 0, 5);
  lprior += normal_lpdf(b_t1 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) | 0, 5));
  lprior += normal_lpdf(b_t2 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t1 | 0, 5));
  lprior += normal_lpdf(b_t3 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t2 | 0, 5));
  lprior += normal_lpdf(b_t4 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t3 | 0, 5));
  lprior += normal_lpdf(b_t5 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t4 | 0, 5));
  lprior += normal_lpdf(b_t6 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t5 | 0, 5));
  lprior += student_t_lpdf(sigma | 3, 0, 11.9)
    - 1 * student_t_lccdf(0 | 3, 0, 11.9);
  lprior += normal_lpdf(sd_1[1] | 0, 1)
    - 1 * normal_lccdf(0 | 0, 1);
  lprior += normal_lpdf(sd_1[2] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[3] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[4] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[5] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[6] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[7] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[8] | 0, 2)
    - 1 * normal_lccdf(0 | 0, 2);
  lprior += normal_lpdf(sd_1[9] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += normal_lpdf(sd_1[10] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += normal_lpdf(sd_1[11] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += normal_lpdf(sd_1[12] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += normal_lpdf(sd_1[13] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += normal_lpdf(sd_1[14] | 0, 5)
    - 1 * normal_lccdf(0 | 0, 5);
  lprior += lkj_corr_cholesky_lpdf(L_1 | 1);
}
8

No. But if you give them distributions based on parameters, then you need the truncation statements. That is, if mu and sigma are parameters, then you need to write:

parameters {
  real mu;
  real<lower=0> sigma;
  vector<lower=LB, upper=UB> [N] alpha;
  ...
model {
  alpha ~ normal(0, 2);  // (1)
  alpha ~ normal(mu, sigma) T[LB, UB];  // (2)
  ...

In case (1), you don’t need the truncation bound because the parameters to the normal are constants. For case (2), you need the truncation because it contributes a non-constant term (that depends on mu and sigma) to the log density.

For example, here you have a parameter b_t4 involved in the lower bound of b_t5.

vector<lower=min(C_1) + b_t4,upper=max(C_1)>[K_t5] b_t5;  // regression coefficients
...
prior += normal_lpdf(b_t5 | 0, 5)
    - 1 * log_diff_exp(normal_lcdf(max(C_1) | 0, 5), normal_lcdf(min(C_1) + b_t4 | 0, 5));

This is correctly doing the adjustment. In general, never multiply by -1, just negate, so that -1 * x becomes -x (it’s faster and uses less memory). You can reduce the lpdf to lupdf to drop unnecessary normalizing constants.

I’m not sure why you’re defining prior as a variable as we have a built in that would let you write this as

b_t5 ~ normal_lpdf(0, 5) T[max(C_1), min(C_1) + b_t4];

You didn’t provide the model, but I take it C_1 is data. I’m generally wary of putting hard boundaries on parameters based on data. For example, if I have a observations y_1, \ldots, y_N, I fit

y ~ normal(mu, sigma);

I wouldn’t declare mu as

real<lower=min(y), upper=max(y)> mu;

It may be that mu values outside the data range are consistent with the data and putting the hard boundaries will give you a different answer even when you expect your posterior mean to lie in that interval.

9

Thank you @Bob_Carpenter

This is very helpful, as always!

Regards