# Offset multiplier initialization

**URL:** <https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712>\
**Category:** General\
**Created:** [February 9, 2021, 7:38pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712 "2021-02-09T19:38:02Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 9, 2021, 7:38pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/1 "2021-02-09T19:38:02Z")

</div>

@jsocolar pointed out that if you initialize an offset-multiplier style non-centered parameters to a really small value you get problems sampling. @jsocolar’s example showed the same problems for me, but when I tried to duplicate it with a manual offset-multiplier implementation, I did not get the problem.

It seems like something not-so-good is going on here.

Here is the offset-multiplier model

```stan
parameters{
  real<lower = 0> sigma;
  real<multiplier=sigma> x;
}
model{
  sigma ~ std_normal();
  x ~ normal(0, sigma);
}

```

Here is the manual offset-multiplier model:

```stan
parameters{
  real<lower = 0> sigma;
  real x_raw;
}
transformed parameters {
  real x = x_raw * sigma;
}
model{
  sigma ~ std_normal();
  x ~ normal(0, sigma);
  target += log(sigma);
}

```

Code to run them is

```
library(tidyverse)
library(cmdstanr)

mod1 = cmdstan_model("mod1.stan")
inits_chain_1 = list(sigma = 1e-20)
fit1 = mod1$sample(chains = 1, init = list(inits_chain_1), iter_sampling = 1000)
fit1$summary()

mod2 = cmdstan_model("mod2.stan")
fit2 = mod2$sample(chains = 1, init = list(inits_chain_1), iter_sampling = 1000)
fit2$summary()

```

You get output like this for the build-in offset-multiplier:

```
> fit1$summary()
# A tibble: 3 x 10
  variable mean median sd mad q5 q95 rhat ess_bulk
  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -9.01e+38 -9.01e+38 0 0 -9.01e+38 -9.01e+38 NA NA
2 sigma 3.12e-20 3.12e-20 0 0 3.12e-20 3.12e-20 NA NA
3 x -1.33e+ 0 -1.33e+ 0 0 0 -1.33e+ 0 -1.33e+ 0 NA NA
# … with 1 more variable: ess_tail <dbl>

```

And the output with a custom offset-multiplier looks like:

```
> fit2$summary()
# A tibble: 4 x 10
  variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -1.57 -1.24 0.987 0.761 -3.64 -0.563 1.00 292. 430.
2 sigma 0.821 0.703 0.594 0.591 0.0667 2.02 1.00 437. 366.
3 x_raw 0.0115 0.0396 0.983 1.01 -1.64 1.59 1.00 523. 660.
4 x 0.00698 0.0104 0.967 0.547 -1.66 1.51 1.00 492. 503.

```

This is pretty repeatable that the custom code doesn’t have a problem with inits but the built in does.

I’m not quite sure where the problem is or what the difference is – posting here if someone knows. @jsocolar said this was happening in a realistic setting (and so figuring out it was a problem probably wasn’t fun or easy and so it’d be nice to fix this if it is fixable or at least know what is the difference)

---

<div class="post-metadata">

**Author:** ![rok\_cesnovar](https://avatars.discourse-cdn.com/v4/letter/r/7bcc69/32.png) [@rok\_cesnovar](https://discourse.mc-stan.org/u/rok_cesnovar)\
**Post date:** [February 9, 2021, 7:47pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/2 "2021-02-09T19:47:12Z")

</div>

Not sure it helps, but the generated C++ for the first case:

```cpp
local_scalar_t__ sigma;
sigma = DUMMY_VAR__;
      
sigma = in__.scalar();
if (jacobian__) {
    sigma = stan::math::lb_constrain(sigma, 0, lp__);
} else {
    sigma = stan::math::lb_constrain(sigma, 0);
}
local_scalar_t__ x;
x = DUMMY_VAR__;

x = in__.scalar();
if (jacobian__) {
    x = stan::math::offset_multiplier_constrain(x, 0, sigma, lp__);
} else {
    x = stan::math::offset_multiplier_constrain(x, 0, sigma);
}
{
    lp_accum__.add(std_normal_lpdf<propto__>(sigma));
    lp_accum__.add(normal_lpdf<propto__>(x, 0, sigma));
}

```

and for the second case:

```cpp
local_scalar_t__ sigma;
sigma = DUMMY_VAR__;

sigma = in__.scalar();
if (jacobian__) {
  sigma = stan::math::lb_constrain(sigma, 0, lp__);
} else {
  sigma = stan::math::lb_constrain(sigma, 0);
}
local_scalar_t__ x_raw;
x_raw = DUMMY_VAR__;

x_raw = in__.scalar();
local_scalar_t__ x;
x = DUMMY_VAR__;

x = (x_raw * sigma);
{
  lp_accum__.add(std_normal_lpdf<propto__>(sigma));
  lp_accum__.add(normal_lpdf<propto__>(x, 0, sigma));
  lp_accum__.add(stan::math::log(sigma));
}

```

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 9, 2021, 7:51pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/3 "2021-02-09T19:51:13Z")

</div>

Thanks and this is the offset multiplier code: [math/offset\_multiplier\_constrain.hpp at develop · stan-dev/math · GitHub](https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/offset_multiplier_constrain.hpp#L80)

Not sure where something is going off the rails. Probably the way to debug this is looking at differences in gradients of those two blocks of C++ tho.

(Edit: to be clear not asking you to check the C++ or anything unless you really want to – the point of the post is to know if anyone else has seen this/knows what is happening before digging in)

---

<div class="post-metadata">

**Author:** ![jsocolar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jsocolar/32/2486_2.png) [@jsocolar](https://discourse.mc-stan.org/u/jsocolar)\
**Post date:** [February 9, 2021, 8:14pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/4 "2021-02-09T20:14:29Z")

</div>

Just to mention that I _think_ the settings where I’ve seen this happen with “real models” involve initializations far from the typical set where gradients are really weak, and so (with nontrivial probability) chains manage to wander down to absurdly low standard deviations during the initial exploration… and (sinister voice) once they go in, they never come out.

I’ve got a model on simulated data where this was happening as frequently as approximately 1 in 4 chains. Even when chains worked, during the first few iterations this model threw a surprising number of warnings about “multiplier is zero but must be positive” and “current metropolis proposal is about to be rejected” (sorry I don’t have the actual warning text in front of me). Switching to the “old” noncentered parameterization squashed all these warnings as well as the sticky-boundary issue.

---

<div class="post-metadata">

**Author:** ![Lu.Zhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lu.zhang/32/8737_2.png) [@Lu.Zhang](https://discourse.mc-stan.org/u/Lu.Zhang)\
**Post date:** [May 9, 2021, 4:16pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/5 "2021-05-09T16:16:34Z")

</div>

I agree with @jsocolar that initialization being far from the high probability region could be the problem. The following plot shows the density plot of the offset-multiplier model. The blue and black region has log-density greater than -15. We can see from the plot that the initials will be far from the high probability region when initializing sigma at a very small value.  
 ![p1](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/9/9403320b9dadd1033c733fe9c75d3bb357a339a0.jpeg)

On the other hand, the high probability region of the manual offset-multiplier has a much better shape. The following plot shows the log-density of the manual offset-multiplier, where the blue and black region has log-density greater than -10.  
 ![p2](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/05ea577f5129fcfd75f5dcdc8445cea30306efc2.jpeg)

The _weird_ shape of the high probability region is probably caused by the singularity of the log-density \log p(x, \sigma) of the offset-multiplier model at \sigma = 0 , where \log p(x, \sigma) \propto -0.5 \sigma^2 - log(\sigma) - 0.5 x^2/\sigma^2. While the manual offset-multiplier model whose log-density follows \log p(x\_{raw}, \sigma) \propto -0.5 \sigma^2 - 0.5 x\_{raw}^2 does not have this problem.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 9, 2021, 5:03pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/6 "2021-05-09T17:03:15Z")

</div>

Oh wait, when I read this issue a couple days ago I thought these were two different models, but now that I look closer I really think these should be exactly the same two log-densities.

This code from the built-in offset-multiplier (jacobian should be true when we’re doing HMC):

```cpp
x = in__.scalar();
if (jacobian__) {
    x = stan::math::offset_multiplier_constrain(x, 0, sigma, lp__);
} else {
    x = stan::math::offset_multiplier_constrain(x, 0, sigma);
}
...
lp_accum__.add(lp__) // this happens later in the code

```

Should be equivalent to:

```cpp
x = (x_raw * sigma);
...
lp_accum__.add(stan::math::log(sigma));

```

So I really think these should be the same log densities? I’m confused. Also a lot of the offset\_multiplier code and whatnot changed in the last few months so maybe some autodiff got improved in develop. I am curious if this would happen there.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 9, 2021, 5:07pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/7 "2021-05-09T17:07:04Z")

</div>

> [@Lu.Zhang](#):
>
> While the manual offset-multiplier model whose log-density follows

I think that density is the one you’d get from a regular non-centered parameterization. Neither of the models in the first post are the regular non-centered parameterization (this is what I missed a couple days ago when I read it).

One of these is offset multiplier and the other should be equivalent to offset multiplier :/.

---

<div class="post-metadata">

**Author:** ![Lu.Zhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lu.zhang/32/8737_2.png) [@Lu.Zhang](https://discourse.mc-stan.org/u/Lu.Zhang)\
**Post date:** [May 9, 2021, 5:37pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/8 "2021-05-09T17:37:26Z")

</div>

Sorry I am a little bit confused. Isn’t it correct that Stan samples in the unconstrained space of (x, sigma) for the offset-multiplier model, and Stan samples in the unconstrained space of (x\_raw, sigma) for the manual offset-multiplier model? The x-axis of the two density plots I gave are different. By the way, I just found that I missed the lower bound of sigma, so I updated the log-density plots in my last post.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 9, 2021, 6:46pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/9 "2021-05-09T18:46:19Z")

</div>

Ooooo, I think you’re right it is initialization. Underlying they are both exposed to the sampler as x\_raw. However, when things are initialized they are initialized on the scale you see, which in the first model is x (not the x\_raw scale). In the second model things are indeed initialized on x\_raw. Not sure that is actually the problem, but could test that by supplying manual inits.

For more on why both models have an x\_raw, what the offset multiplier thing does behind the scenes in sample on an x\_raw to avoid that by doing a change of variables transformation from an x\_raw to x that should mimic what the second model does.

That’s mostly baked into:

```
x = stan::math::offset_multiplier_constrain(x, 0, sigma, lp__);

```

There is the x on the input but that is different than the x on the output (the input should be equivalent to x\_raw).

---

<div class="post-metadata">

**Author:** ![Lu.Zhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lu.zhang/32/8737_2.png) [@Lu.Zhang](https://discourse.mc-stan.org/u/Lu.Zhang)\
**Post date:** [May 9, 2021, 6:55pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/10 "2021-05-09T18:55:07Z")

</div>

> [@bbbales2](#):
>
> For more on why both models have an x\_raw, what the offset multiplier thing does behind the scenes in sample on an x\_raw to avoid that by doing a change of variables transformation from an x\_raw to x that should mimic what the second model does.
> 
> That’s mostly baked into:
> 
> ```no-highlight
> x = stan::math::offset_multiplier_constrain(x, 0, sigma, lp__);
> 
> ```

I see, that is interesting! Thank you so much!

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 9, 2021, 6:59pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/11 "2021-05-09T18:59:06Z")

</div>

> [@bbbales2](#):
>
> when things are initialized they are initialized on the scale you see

Wait, I say that with such confidence but it isn’t true lol. Nevermind whoops.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [May 9, 2021, 7:00pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/12 "2021-05-09T19:00:02Z")

</div>

(Randomly initialized things are initialized on the unconstrained scale)

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [May 10, 2021, 7:21am UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/13 "2021-05-10T07:21:36Z")

</div>

Wow, that looks weird.

If that helps anybody, I tried building a Stan model that should match all of the operations and their order as in the offset-multiplier (mod1) and it still does not match the (bad) results in mod1:

`mod3.stan`:

```stan
parameters{
  real<lower = 0> sigma;
  real x_raw;
}
transformed parameters {
}
model{
  real x;
  target += log(sigma);
  x = fma(sigma, x_raw, 0);
  sigma ~ std_normal();
  x ~ normal(0, sigma);
}

```

So the only remaining difference IMHO is that in mod1 the constraining transform overwrites the `x` local variable while in mod3 we create a new local variable to hold the result… Which should not matter much?

BTW, I’ve always found the offset-multiplier way of doing stuff confusing as a) it requires me to use `sigma` at two potentially distant places in the code, increasing risk I will change one but not the other in the future and b) I always have to think hard whether the multiplier is from unconstrained-to-constrained or the other way around. In fact, our own manual used to have the wrong way! (contrast [2.26 manual](https://mc-stan.org/docs/2_26/reference-manual/affinely-transformed-scalar.html) with the [current version](https://mc-stan.org/docs/reference-manual/transforms.html#affinely-transformed-scalar), see also [docs issue #356](https://github.com/stan-dev/docs/issues/356) )

In addition the offset-multiplier logic (both manual and controled by Stan) adds unnecessary operations to the model: `log(sigma)` is added in the Jacobian correction, only to be subtracted in `x ~ normal(0, sigma)`, with a penalty for both performance and numerical precision, so I think the best way to code non-centering is the IMHO simpler:

```stan
parameters{
  real<lower = 0> sigma;
  real x_raw;
}
transformed parameters {
  real x = x_raw * sigma;
}
model{
  sigma ~ std_normal();
  x_raw ~ std_normal();
}

```

---

<div class="post-metadata">

**Author:** ![Lu.Zhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/lu.zhang/32/8737_2.png) [@Lu.Zhang](https://discourse.mc-stan.org/u/Lu.Zhang)\
**Post date:** [May 10, 2021, 7:46pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/14 "2021-05-10T19:46:00Z")

</div>

Hmm~ I still suspect that the way of generating initials of the two models are different.  
If you check the first several samples of x and \sigma for the two models

```
mod1 = cmdstan_model("mod1.stan")
inits_chain_1 = list(sigma = 1e-7)
fit1 = mod1$sample(chains = 1, init = list(inits_chain_1), 
               iter_sampling = 1000, save_warmup = TRUE)
draws <- fit1$draws(inc_warmup = TRUE)
rbind(draws[1:6,,"x"], draws[1:6,,"sigma"])

mod2 = cmdstan_model("mod2.stan")
fit2 = mod2$sample(chains = 1, init = list(inits_chain_1), 
                   iter_sampling = 1000, save_warmup = TRUE)
draws2 <- fit2$draws(inc_warmup = TRUE)
rbind(draws2[1:6,,"x"], draws2[1:6,,"sigma"])

```

you will get output like this:

```
> rbind(draws[1:6,,"x"], draws[1:6,,"sigma"])
          [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.2110230 0.2110230 -1.30120e-02 -7.56263e-03 -1.14654e-02 -1.61252e-02
[2,] 0.0000001 0.0000001 1.38811e-07 8.30654e-08 1.32924e-07 2.14542e-07
> rbind(draws2[1:6,,"x"], draws2[1:6,,"sigma"])
            [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 2.27659e-08 -8.2056e-08 6.52913e-06 -6.68192e-06 0.000232000 0.000181418
[2,] 1.00000e-07 4.5713e-07 5.38727e-06 5.36030e-06 0.000327161 0.003488220

```

If both models initialize on the unconstrained space of x\_{raw} and \sigma, why the scale of x are so different for the two models.

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [May 10, 2021, 8:12pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/15 "2021-05-10T20:12:31Z")

</div>

In case anyone wants to get confused a bit more, here’s a traceplot of the first 200 warmup iterations for a single chain and `mod1` and `mod2` with a fixed seed but unspecified inits:

 ![trace](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/72e3a3c31752aaa3cac95549a72fa7c7de8fae36.png)

Edit with further info:

- Any common seed without inits yields a picture as above
- Specifying only `sigma` or only `x` and `x_raw` as inits yields different initial points and initial traces (as for @Lu.Zhang)
- Specifying `sigma` and `x==x_raw/sigma` again yields the same inits and slowly diverging traces.

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [May 10, 2021, 8:40pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/16 "2021-05-10T20:40:18Z")

</div>

Also, @martinmodrak wins a prize:

The other formulations start to do very weird things around `sigma==1e-10`:

 ![trace_1e-10](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/4/4308e18520e5400b9a49cb10584b7e8d05bdf4c4.png)

Another edit with further info:

`unit` corresponds to

> [@martinmodrak](#):
>
> ```stan
> parameters{
> real<lower = 0> sigma;
> real x_raw;
> }
> transformed parameters {
> real x = x_raw * sigma;
> }
> model{
> sigma ~ std_normal();
> x_raw ~ std_normal();
> }
> 
> ```

Final edit:

left vs right: initialize only sigma vs initialize sigma and x/x\_raw  
top to bottom: sigma from 1e-20 to 1e-4  
colors: models  
lines: traceplots in the sigma,|x| plane

black lines: common initial values (if applicable)

 ![trace](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0af8f704c478a50591ce3c86472c10c87f346f1d.png)

Summary:

The manual version only appeared better because the initial points were not the same if only sigma was specified. Actually, it’s not really better than the offset-multiplier version. The reason appears to be numerical, as @martinmodrak’s version works fine for all (tested) values of sigma.

---

<div class="post-metadata">

**Author:** ![Dalton](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dalton/32/4287_2.png) [@Dalton](https://discourse.mc-stan.org/u/Dalton)\
**Post date:** [April 15, 2025, 11:17pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/17 "2025-04-15T23:17:13Z")

</div>

Coming in 4 years later to chime in that I recently ran into this issue in wild when attempting to fit a logit-normal hierarchical binomial model even with relatively a zero-avoiding prior for sigma. I’m able to achieve convergence by either setting a very small non-zero lower bound on sigma (1e-10) or by switching to a manual non-centered form. But with a lower bound of zero and the affine transform, this model pretty reliably fails.

Here is a reproducible and realistic example of this behavior:

```no-highlight
data {
  int<lower = 0> I;
  array[I] int N, y;
  real<lower = 0> sigma_lower_bound;
}

parameters {
  real mu;
  real<lower = sigma_lower_bound> sigma;
  vector<offset=mu, multiplier=sigma>[I] logit_p;
  // vector[I] logit_p_raw;
}

transformed parameters{
  // vector[I] logit_p = mu + logit_p_raw * sigma;
}

model{
  mu ~ student_t(3, 0, 2.5);
  sigma ~ loglogistic(1, 4);
  logit_p ~ normal(mu, sigma);
  // logit_p_raw ~ std_normal();
  
  y ~ binomial_logit(N, logit_p);
}

```

In R:

```no-highlight
library(cmdstanr)
library(bayesplot)

I = 1000
p = plogis(rnorm(I, mean = 0.69, sd = 0.5))
N = rpois(I, rgamma(I, 20, 1))
y = rbinom(I, N, p)

binom_hier <- cmdstan_model("binomial_hier.stan")

binom_hier_fit <- binom_hier$sample(data = list(I = I, 
                                                N = N, 
                                                y = y, 
                                                sigma_lower_bound = 0),
                                  chains = 4,
                                  parallel_chains = 4,
                                  refresh = 100,
                                  iter_warmup = 500,
                                  iter_sampling = 500,
                                  save_warmup = T)

```

Here is the result from running this model:

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

In this example, the initial values weren’t all that close to the boundary. Chain 1 was 1.3, Chain 2 0.5, Chain 3 0.31 and Chain 4: 1.4. By iteration 5 Chain 2&3 go down to 1e-16 and never recover, but Chain 4 briefly flirts with 1e-10 before bouncing back. It’s like an event horizon down there somewhere.

---

<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:** [April 16, 2025, 10:07pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/18 "2025-04-16T22:07:34Z")

</div>

> [@bbbales2](#):
>
> However, when things are initialized they are initialized on the scale you see

Other than rejecting for 0 values, you should be able to exactly duplicate behavior between the two approaches if you use your own custom initialization.

The declaration

```stan
real<offset=mu, multiplier=sigma> alpha;

```

means that the “constrained” variable `alpha` is transformed to the “unconstrained” value `(alpha - mu) / sigma`. This means that Stan’s default uniform random initialization, which is generically \text{uniform}(-2, 2) in unconstrained space, becomes uniform on (\mu - 2 \sigma, \mu + 2\sigma) in the constrained space where the model is defined. If you manually initialize this way, the behavior should match.

> [@Dalton](#):
>
> Here is a reproducible and realistic example

This is just one model, so I assume you want us to edit it in the obvious way for the comparison (a reproducible example would require two separate models so they could be compared).

The two models are not equivalent because of the initialization.

P.S. If you have a problem where a lower bound on a scale like `sigma` feels like a solution, it’s signaling a bigger problem with your model that the lower bound is just sweeping under the carpet (almost always that the data’s consistent with zero variance). If you really want to keep things away from zero, the loglogistic should do that on its own as the density goes to zero as the value goes to zero.

Unless you really think your `mu` values could be relatively large, you don’t need a wide-tailed distribution like Student-t with 3 degrees of freedom in the prior. Your 99% intervals are (-15, 15) and the 99.9% are (-30, 30), but the tail goes out into values like 100 not being unreasonable. If -30 isn’t a reasonable value, that prior’s too wide. This can wind up causing a problem with fits if you don’t have much data. If you have a ton of data, it should be fine.

---

<div class="post-metadata">

**Author:** ![Dalton](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dalton/32/4287_2.png) [@Dalton](https://discourse.mc-stan.org/u/Dalton)\
**Post date:** [April 16, 2025, 11:54pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/19 "2025-04-16T23:54:56Z")

</div>

> P.S. If you have a problem where a lower bound on a scale like `sigma` feels like a solution, it’s signaling a bigger problem with your model that the lower bound is just sweeping under the carpet (almost always that the data’s consistent with zero variance). If you really want to keep things away from zero, the loglogistic should do that on its own as the density goes to zero as the value goes to zero.

But the data aren’t consistent with zero variance. It was simulated with non-zero variance. I just reran the model using I = 500, and 8 chains. One got stuck (another flirted with disaster in warm-up and recovered). The `lp__` for the chain that got stuck at zero is a lot lower than for the other 7 chains.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/8/a/8aff5833846ece99266757954aa0faa4d1eb79f4.png)  
Using `glm()` in R on that same data tell me that the dispersion parameter is 2.

> Other than rejecting for 0 values, you should be able to exactly duplicate behavior between the two approaches if you use your own custom initialization.

I’m not convinced this is the initial values, because in the example I posted yesterday, the first iteration values for `sigma` were all reasonably far from zero. Now I know these aren’t the randomly generated initial values however. So I’m currently working on finding a set that will reproduce this behavior exactly. From trial and error, I can say with the logit-normal hierarchical model, it seems like I’m more likely to run into this issue when I increase the number of draws from the hyperparameter (that is increase I). So maybe the issue is one of the unconstrained `logit_phi[i]` values being initialized in a way that doesn’t play nicely with the unconstrained `sigma` and `mu` values? I’ll see if I can play around and find some initial values that reliably reproduce this.

---

<div class="post-metadata">

**Author:** ![Dalton](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dalton/32/4287_2.png) [@Dalton](https://discourse.mc-stan.org/u/Dalton)\
**Post date:** [April 17, 2025, 9:19pm UTC](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712/20 "2025-04-17T21:19:29Z")

</div>

Okay, so reporting back in with perhaps stronger evidence of something funky happening with the affine transform in this logit-normal hierarchical model. And while it does seem to have something to do with initial values, it doesn’t appear to be anything in particular about low initial values for `sigma`.

As suggested I am comparing two “identical” models fit to the same data with the same initial values and the same seed.

As @Bob_Carpenter suggested here are two separate scripts:  
`binomial_hier_affine.stan`

```no-highlight
data {
  int<lower = 0> I;
  array[I] int N, y;
}
parameters {
  real mu;
  real<lower = 0> sigma;
  vector<offset=mu, multiplier=sigma>[I] logit_p;
}
model{
  mu ~ student_t(3, 0, 2.5);
  sigma ~ loglogistic(1, 4);
  logit_p ~ normal(mu, sigma);

  y ~ binomial_logit(N, logit_p);
}

```

and  
`binomial_hier_raw.stan`

```no-highlight
data {
  int<lower = 0> I;
  array[I] int N, y;
}
parameters {
  real mu;
  real<lower = 0> sigma;
  vector[I] logit_p_raw;
}
transformed parameters{
  vector[I] logit_p = mu + logit_p_raw * sigma;
}
model{
  mu ~ student_t(3, 0, 2.5);
  sigma ~ loglogistic(1, 4);
  logit_p_raw ~ std_normal();
  
  y ~ binomial_logit(N, logit_p);
}

```

In R we simulate overdispersed binomial data like so:

```no-highlight
I = 500
p = plogis(rnorm(I, mean = 0.69, sd = 0.5))
N = rpois(I, rgamma(I, 20, 1))
y = rbinom(I, N, p)
x = rbinom(I, N, 0.69)

summary(glm(cbind(y, N-y) ~ NULL, family = "quasibinomial"))

```

We can see that the data are consistent with extrabinomial variation by checking with `glm()`

```no-highlight
(Dispersion parameter for quasibinomial family taken to be 1.935123)

```

Here is function to generate initial values ensuring that the both models receive the same initial values for both `sigma` and `logit_p`/`logit_p_raw` parameters

```no-highlight
init_binom <- function(mu_in, logit_p_raw_in, sigma_in, rand = T){
  inits = list()
  for (i in 1:length(sigma_in)){
    if(rand){
      logit_p_raw = runif(length(logit_p_raw_in[,1]), -2, 2)
      sigma_in = exp(runif(1, -2, 2))
      mu_in = runif(1, -2, 2)
      inits[[i]] = list(mu = mu_in, sigma = sigma_in, 
                        logit_p = mu_in + logit_p_raw*sigma_in,
                        logit_p_raw = logit_p_raw)
    } else{
      logit_p_raw = logit_p_raw_in
      inits[[i]] = list(mu = mu_in[i], sigma = sigma_in[i], 
                        logit_p = mu_in[i] + logit_p_raw[,i]*sigma_in[i],
                        logit_p_raw = logit_p_raw[,i])
    }
  }
  return(inits)
}

```

And some initial values tests, one deliberately checking low initial values of `sigma` and holding other parameters fixed and the other checking with random values that I think should be similar to how Stan generates them.

```no-highlight
mu_init = rep(0.69, 8)
sigma_chain = c(1e-21, 1e-15, 1e-10, 1e-6, 1e-3, 0.01, 0.1, 1)
logit_p_raw_init = matrix(nrow = I, ncol = 8)
for (i in 1:8)
  logit_p_raw_init[,i] = rep((0.7 - mu_init[i]) / sigma_chain[i], I)

set.seed(125)

inits_nonrandom = init_binom(mu_init, logit_p_raw_init, sigma_chain, F)
inits = init_binom(mu_init, logit_p_raw_init, sigma_chain, T)

```

Running the models with non-random inits shows that the affine transform reliably fails when `sigma` is initialized to low values. The manual version also struggled (with one chain giving me treedepth warnings) but all of the chains converged.

First a check that the models are initializing to the same values as suggested earlier in this thread:

```no-highlight
init_test_affine = binom_hier_fit_affine$draws(variables = c("mu", "sigma", "logit_p[1]", "lp__"), inc_warmup = T)
init_test_raw = binom_hier_fit_raw$draws(variables = c("mu", "sigma", "logit_p[1]", "lp__"), inc_warmup = T)

```

```no-highlight
> init_test_affine[1:6,,"sigma"]
# A draws_array: 6 iterations, 8 chains, and 1 variables
, , variable = sigma

         chain
iteration 1 2 3 4 5 6 7 8
        1 1e-21 4.7e-90 4.4e-159 1.0e-06 0.1210 0.26 0.34 0.80
        2 1e-21 4.7e-90 4.4e-159 1.0e-06 0.1210 0.26 0.34 0.80
        3 1e-21 3.6e-67 4.4e-159 6.7e-04 0.1210 0.26 0.34 0.80
        4 1e-21 1.3e-63 4.7e-159 2.3e-264 0.1210 0.26 0.34 0.80
        5 1e-21 1.3e-63 6.6e-159 2.2e-264 0.0026 0.38 0.45 0.74
> init_test_raw[1:5,,"sigma"]
# A draws_array: 5 iterations, 8 chains, and 1 variables
, , variable = sigma

         chain
iteration 1 2 3 4 5 6 7 8
        1 1.0e-21 2.8e-13 1.0e-10 1.0e-06 0.1210 0.26 0.34 0.80
        2 1.0e-21 2.8e-13 1.0e-10 1.0e-06 0.1210 0.26 0.34 0.80
        3 1.5e-21 2.8e-13 9.4e-09 6.7e-04 0.1210 0.26 0.34 0.80
        4 1.7e-20 9.5e-242 8.5e-153 2.3e-264 0.1210 0.26 0.34 0.80
        5 2.5e-45 1.2e-241 1.1e-152 3.1e-264 0.0026 0.38 0.45 0.74
> init_test_affine[1:5,,"logit_p[1]"]
# A draws_array: 5 iterations, 8 chains, and 1 variables
, , variable = logit_p[1]

         chain
iteration 1 2 3 4 5 6 7 8
        1 0.71 31233 4097 0.7 0.44 0.64 0.55 1.3
        2 0.71 31233 4097 0.7 0.44 0.64 0.55 1.3
        3 0.71 -9704 4097 330.3 0.44 0.64 0.55 1.3
        4 0.71 3404 3699 225.2 0.44 0.64 0.55 1.3
        5 0.71 3404 3524 -18.9 0.64 0.97 1.04 1.7
> init_test_raw[1:5,,"logit_p[1]"]
# A draws_array: 5 iterations, 8 chains, and 1 variables
, , variable = logit_p[1]

         chain
iteration 1 2 3 4 5 6 7 8
        1 0.7 -322 0.7 0.7 0.44 0.64 0.55 1.3
        2 0.7 -322 0.7 0.7 0.44 0.64 0.55 1.3
        3 -3.7 -322 325.2 330.3 0.44 0.64 0.55 1.3
        4 176.4 -133 125.5 225.2 0.44 0.64 0.55 1.3
        5 492.3 209 -56.6 43.0 0.64 0.97 1.04 1.7

```

That’s a qualified yes, though I’m not sure what’s going on with chains 2 and 3.

Traceplots. For the affine transform model, all four chains with `sigma` initialized \< 1e-5 got stuck with `sigma` near-zero, but oddly they all bottomed out far below the initialized values before somewhat recovering.

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

Here are the traces for manual non-centering model:

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

Now for something very odd with the the random inits. I was able to pretty quickly hit upon a combination that resulted in a single chain getting stuck at zero for sigma using the affine transform version:

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

But the initial value for chain 4 was pretty high for `sigma`. The value provided was 2.26, which got turned into 0.84 for the first few warm-up iterations:

```no-highlight
> unlist(c(inits[[1]][2], inits[[2]][2], inits[[3]][2], inits[[4]][2],
+ inits[[5]][2], inits[[6]][2], inits[[7]][2], inits[[8]][2]))
    sigma sigma sigma sigma sigma sigma sigma sigma 
0.2335458 6.3053346 4.5656075 2.2926901 0.2683244 0.4980496 0.2839859 4.6789906 
> init_test_affine[1:5,,"sigma"]
# A draws_array: 5 iterations, 8 chains, and 1 variables
, , variable = sigma
         chain
iteration 1 2 3 4 5 6 7 8
        1 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        2 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        3 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        4 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        5 8.0e-08 0.38 0.148 0.39 8.8e-02 0.29 0.00021 1.5e-15
> init_test_raw[1:5,,"sigma"]
# A draws_array: 5 iterations, 8 chains, and 1 variables
, , variable = sigma

         chain
iteration 1 2 3 4 5 6 7 8
        1 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        2 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        3 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        4 2.3e-01 2.42 0.026 0.84 2.3e-13 0.43 0.14313 2.8e-01
        5 8.0e-08 0.38 0.148 0.39 8.8e-02 0.29 0.00021 4.0e-15

```

Plotting the paths of the two models for the first 20 iterations (iteration number labelled) and you can see them track each other pretty well at first. But then chain 4 gets stuck for the affine transform. The other two chains (5 and 8) are displayed because for these chains values also started to get down near zero for `sigma` but then recovered (none of the other chains ever got quite so low). So I’m really not clear what combination of conditions is resulting in the affine transform performing so poorly with this model.

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

In summary, there is something different about the affine transformation compared to the manual non-centered transform that results in different behavior. Identical initial values and seeds will diverge and it’s not solely a feature of poor initial conditions in terms of `sigma` being set to zero. However, what’s going on here is a sort of triple transformation (lower bound for `sigma`, affine/non-centering for the `logit_p`’s, and then `p` for binomial. So maybe something pathological in those combinations of constraints?

[Next page](https://discourse.mc-stan.org/t/offset-multiplier-initialization/20712.md?page=2)