# Avoid extreme values or putting boundaries

**URL:** <https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687>\
**Category:** Modeling\
**Created:** [September 21, 2024, 10:40pm UTC](https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687 "2024-09-21T22:40:27Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![dmp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dmp/32/16409_2.png) [@dmp](https://discourse.mc-stan.org/u/dmp)\
**Post date:** [September 21, 2024, 10:40pm UTC](https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687/1 "2024-09-21T22:40:27Z")

</div>

hi all, I am using stan to build a Bayesian hierarchical model where individual-level parameters are a function of their individual-characteristics. The structure of my model is similar to the following descriptions. I am trying to use cholesky\_factor\_corr() to reparameterize it. My question is, how should I put constraints to avoid my individual-level parameters gamma\_i to have extremely big/small values? I have already put lower bounds for the scale e.g., lambda\_raw, etc., however, it still produces values as small as -10^100. Should I change priors ?

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/e/f/efe0fabc03ecdeea8b01801af9da638b10640898.png)

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/2/a/2ad90e0c91bc6ce899627a859d42665a35562eb0.png)

```no-highlight
data{ //with reparameterization.
    //other irrelevant part...

    matrix[ncust,nchar] Demo;
}

transformed data {
   matrix[nchar+1,ncust] Demo_tran;
   row_vector[ncust] ones = rep_row_vector(1,ncust);
   Demo_tran = append_row(ones,transpose(Demo));
}

parameters{
    vector<lower=-1,upper=1>[3*(nchar+1)] lambda_raw;
    vector<lower=0,upper=1>[3*(nchar+1)] sigma_lambda;
    cholesky_factor_corr[3*(nchar+1)] L_lambda;

    matrix<lower=-1,upper=1>[3,ncust] gamma_raw;

    vector<lower=0,upper=3>[3] tau_gamma;
    cholesky_factor_corr[3] L_gamma;
    //...other irrelevant part
}

transformed parameters {
    vector[3*(nchar+1)] lambda_vec;
    lambda_vec = sigma_lambda.*(L_lambda*lambda_raw);

    matrix[3,nchar+1] Lambda_mat;
    Lambda_mat = to_matrix(lambda_vec,3,nchar+1);

    matrix[3,3] L_V_gamma = diag_pre_multiply(tau_gamma,L_gamma);
    matrix[3,ncust] mu_gamma;
    mu_gamma = Lambda_mat * Demo_tran;

    matrix<lower=-700,upper=50>[3,ncust] Gamma = mu_gamma + L_V_gamma*gamma_raw;
   
//.... other irrelevant part

}

model{
    lambda_raw ~ std_normal();
    sigma_lambda ~ cauchy(0,5);
    L_lambda ~ lkj_corr_cholesky(2);

    tau_gamma ~ cauchy(0,2.5);
    L_gamma ~ lkj_corr_cholesky(2);

    to_vector(gamma_raw) ~ std_normal();

///....
    }

}

```

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 27, 2024, 8:37pm UTC](https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687/2 "2024-09-27T20:37:46Z")

</div>

Sorry this didn’t get answered sooner—I just saw it as it got promoted to the unanswered category.

I would suggest removing those hard upper and lower interval constraints on parameters like `lambda_raw` and `sigma_lambda`—if probability bunches up against the boundary, it will skew parameter estimates vs. softer parameterizations.

Similarly, if you’re getting estimates of parameters so near zero that they’re close to underflowing, it means your model is consistent with the value being zero. You might want to just set such parameters to zero if that’s where most of the mass is. Otherwise, if you really do want values of `1e-100`, then it might help to reparameterize so that the unconstrained values are more unit scaled.

The following constraint,

```stan
matrix<lower=-700,upper=50>[3,ncust] Gamma = mu_gamma + L_V_gamma*gamma_raw;

```

is probably not doing what you want it to do. This will evaluate whether the right hand side satisfies the bound, and if not, just reject the draw. This can reduce HMC to a form of super inefficient rejection sampling. Instead, you need to jointly constraint the variables, for example, setting the bounds on `mu_gamma` to be `<lower=-700 - L_V_gamma*gamma_raw, upper = 50 - L_V_gamma*gamma_raw>` will ensure `Gamma` follows the right constraints without any rejection. But again, I’d probably try to avoid upper and lower-bounding this altogether.

This, `sigma_lambda ~ cauchy(0,5);`, is a very wide prior that does pretty much nothing given your constraint, `vector<lower=0,upper=1>[3*(nchar+1)] sigma_lambda;`. You might as well just let that be uniform. A more stable alternative would be to use a weakly informative half-normal prior with only a `lower=0` constraint on the parameter.

To figure out what’s going on with parameters on the scale of `1e-100` after warmup, you want to figure out how they’re getting there. Is it that the `lambda_vec` values want to be near zero and you get there by pushing `lambda_raw` to zero? When the fit crashes, it’s an opportunity to investigate and sharpen up modeling assumptions. Samplers work much better when the model’s well specified for the data—when the model and data don’t match well, the sampler has to struggle with challenging geometry to get an answer.

P.S. You can put these on a single line and reorg to get rid of need for parens (arithmetic is left associative in Stan, so `a * b * c` is evaluated as `(a * b) * c`:

```stan
vector[3*(nchar+1)] lambda_vec = L_lambda * lambda_raw .* sigma_lambda;

```

---

<div class="post-metadata">

**Author:** ![dmp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dmp/32/16409_2.png) [@dmp](https://discourse.mc-stan.org/u/dmp)\
**Post date:** [September 28, 2024, 7:18pm UTC](https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687/3 "2024-09-28T19:18:02Z")

</div>

Hi, this is actually very helpful. I didn’t know that we could do a joint constraint on the variables. and I also vaguely remember I read somewhere that uniform can be a improper prior (I forget the source). I will definitely try on this.  
I ran into another issue when I was estimating the model. So in the paper I am following, there are many levels of hierarchy as shown in below:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/6/2/62cd791323de03517977e98760bec59db48cbb85.png)

where the Q\_ijt, for t =0,1,…, t, and Q\_bar\_ijt, v^2\_ijt can be thought of transformed parameters based on the data and the parameters Q\_ij, and sigma^2\_D\_i.

here is the function that specifies the relationship (for now we can ignore the term that contains the \sigma^2\_F, and just focus on the sigma^2\_D\_I for illustration):

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/7/2/724fed2354f3a43a99333787a1951302af34f486.png)

I have also tried to reparameterize the Q\_ijt, but it seems that the draws are still being inefficient. The model was fine when I only add heterogeneity for Q\_j. The minute I try to add heterogeneity in other parameters (that could interact with the variance), I ran into issues.

Here is the part of the code that is relevant, where M\_cus corresponds to the Q\_ijt,  
and this line (in the transformed parameter block) corresponding to the equation (15)

```no-highlight
 Ms[up] = fmin(fmax(Ms_bar[up] + Ms_std[up]*sqrt_Vbar[up],-1e6),1e6);

```

that these two lines (in the transformed parameter block) corresponding to equation (13) and (14)

```no-highlight
       Vs_bar[up] = (K[up]*square(v0_safe)*V)/square(denominator);
        sqrt_Vbar[up] = sqrt(Vs_bar[up]+1e-10);
...
        Ms_bar[up] = exp(log_vs-log_v0)*m0 + exp(log_vs-log(V)+log(K[up]))*M;

```

I have used the safe\_v0 to avoid the nan instances.

```no-highlight
data{
///...
    matrix[ncust,ndemo] Demo;
}

parameters {
    vector[ncust] M_cus_std;
    real M_mean;M_mean;
    real<lower=0> M_sig;

    vector[ndemo] m_par;

    real<lower=0> r;
    vector[ndemo] r_par;
    array[nupdt] real Ms_std;
    real<lower=0> V;
//...    
}

transformed parameters {
    vector[ncust] M_cus = M_mean + M_sig * M_cus_std + Demo * m_par;
    vector[ncust] r_cus = r + Demo * r_par;

    array[nupdt] real Ms;
    array[nupdt] real Ms_bar;

    array[nupdt] real<lower= 0> Vs_bar;
    array[nupdt] real <lower= 0> Vs;
    array[nupdt] real <lower= 0> sqrt_Vbar;
    
    for (up in 1:nupdt){
        int cus = Uid_updt[up];
        real M = M_cus[cus];

        real m0 = (Itd_updt[up]==1)? M_mean:Ms[up-1];
        real v0 = (Itd_updt[up]==1)? ini_var:Vs[up-1];
        
        real v0_safe = fmax(v0,1e-8); //Ensure v0 is not too close to zero
        real denominator = fmax(v0_safe*K[up]+V,1e-8); //prevent division from very small numbers

        Vs[up] = (v0_safe*V)/denominator;
        Vs_bar[up] = (K[up]*square(v0_safe)*V)/square(denominator);

        sqrt_Vbar[up] = sqrt(Vs_bar[up]+1e-10);
        real log_vs = log(Vs[up]);
        real log_v0 = log(v0_safe);
        Ms_bar[up] = exp(log_vs-log_v0)*m0 + exp(log_vs-log(V)+log(K[up]))*M;

        Ms[up] = fmin(fmax(Ms_bar[up] + Ms_std[up]*sqrt_Vbar[up],-1e6),1e6);
   }
}

model{
    
    M_mean ~ normal(0, 1);
    M_sig ~ normal(0,2.5);//exponential(3);
    M_cus_std ~ std_normal();//normal(0, 1);

    m_par ~ normal(1,1);

    r ~ lognormal(log(1.5),0.5);
    r_par ~ normal(1,1);

    V ~ lognormal(log(4),0.5);
    ini_var ~ lognormal(log(10),1);

    Ms_std ~ normal(0,1);
    intcept ~ normal(2,1);
    restpar ~ normal(0,0.5);

```

I got the following error message, and the R\_hat is being really big.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/5/d/5d2c10cc9543603d8aa6ef4e37f068ada57d1359.jpeg)

I am getting desparate after trying many different ways and it still didn’t work. Any advices are deeply appreciated!!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [October 3, 2024, 9:53pm UTC](https://discourse.mc-stan.org/t/avoid-extreme-values-or-putting-boundaries/36687/4 "2024-10-03T21:53:42Z")

</div>

> [@dmp](#):
>
> I read somewhere that uniform can be a improper prior

Uniform is only proper on a bounded interval. Otherwise the integral diverges.

> [@dmp](#):
>
> The model was fine when I only add heterogeneity for Q\_j. The minute I try to add heterogeneity in other parameters (that could interact with the variance), I ran into issues.

This can potentially be the model and not the way you coded it. When the are a lot of things to estimate you need a lot of data.

Doing hard mins and maxes like you have is a bad idea in Stan because it disconnects the gradients whenever the condition returns the constant `-1e6` or `1e6`. Models usually fit better with softer constraints and usually make more sense unless there’s a hard physical constraint. Whenever you have to start fiddling with epsilons, there’s usually a better solution.

The divergences usually come from varying curvature in hierarchical models when the data is consistent with complete pooling (that drives the hierarchical variance parameter toward zero, which requires very small step sizes, leading to most of the divergences).

> [@dmp](#):
>
> I am getting desparate after trying many different ways and it still didn’t work.

Then you might be ready to try simulating data from your model and fitting it. If that doesn’t work, you know something’s wrong other than the model being misspecified. It helps to start from a simpler model and build up, but it looks like you’re already doing that.

If you look at trace plots of those values with high R-hat values, are they separated? If so, it may be a multimodal posterior.