# Interval-censored normal model: integration vs estimation

**URL:** <https://discourse.mc-stan.org/t/interval-censored-normal-model-integration-vs-estimation/39817>\
**Category:** Modeling\
**Created:** [July 1, 2025, 6:04am UTC](https://discourse.mc-stan.org/t/interval-censored-normal-model-integration-vs-estimation/39817 "2025-07-01T06:04:52Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mdonoghoe](https://avatars.discourse-cdn.com/v4/letter/m/8dc957/32.png) [@mdonoghoe](https://discourse.mc-stan.org/u/mdonoghoe)\
**Post date:** [July 1, 2025, 6:04am UTC](https://discourse.mc-stan.org/t/interval-censored-normal-model-integration-vs-estimation/39817/1 "2025-07-01T06:04:52Z")

</div>

Suppose y\_i \sim N(\mu, \sigma) (i = 1, \ldots, N) but we observe y^\*\_i = \lfloor y\_i \rfloor.

A model for inference about \mu, \sigma can be based on integration, i.e.:

p(y^\*\_i \mid \mu, \sigma) = \int^{y^\*\_i + 1}\_{y^\*\_i} \mathrm{normal}(y \mid \mu, \sigma) dy

In Stan code:

```stan
data {
  
  int<lower=1> N;
  array[N] int y_star;
  
}

parameters {
  
  real mu;
  real<lower=0> sigma;
  
}

model {
  
  for (n in 1:N) {
    
    target += log_diff_exp(normal_lcdf(y_star[n] + 1 | mu, sigma), 
                           normal_lcdf(y_star[n] | mu, sigma));
    
  }
  
}

```

However, because I want to eventually adapt this idea to a different problem, I would like to “estimate” the unobserved y\_i, similar to the approach taken in the [Stan user’s guide](https://mc-stan.org/docs/stan-users-guide/truncation-censoring.html#estimating-censored-values).

In that example there is a fixed censoring value, so it uses a simple constraint in the parameter declaration. I note that [it is possible to set varying constraints](https://mc-stan.org/docs/stan-users-guide/reparameterization.html#varying-upper-and-lower-bounds), but again because of my planned extension to another case, I would like to do it in a more complicated way: basically using an inverse transform sample:

y\_i = \mu + \sigma \Phi^{-1}\left( \Phi(\frac{y^\*\_i - \mu}{\sigma}) + u\_i (\Phi(\frac{y^\*\_i + 1 - \mu}{\sigma}) - \Phi(\frac{y^\*\_i - \mu}{\sigma})) \right) where u\_i \sim U(0,1).

If I want y\_i \sim N(\mu, \sigma), I need to compute the Jacobian:

\frac{dy\_i}{du\_i} = \sigma \frac{\Phi\left(\frac{y^\*\_i + 1 - \mu}{\sigma}\right) - \Phi\left(\frac{y^\*\_i - \mu}{\sigma}\right)}{\phi\left(\frac{y\_i - \mu}{\sigma}\right)}

The other parameters aren’t transformed and the Jacobian is triangular, so its log-determinant is just the log of the value above. In Stan:

```stan
data {
  
  int<lower=1> N;
  array[N] int y_star;
  
}

parameters {
  
  real mu;
  real<lower=0> sigma;
  array[N] real<lower=0,upper=1> u;
  
}

transformed parameters {
  
  array[N] real y;
  
  for (n in 1:N) {
    
    y[n] = mu + sigma * inv_Phi(normal_cdf(y_star[n] | mu, sigma) + 
      u[n] * (normal_cdf(y_star[n] + 1 | mu, sigma) - normal_cdf(y_star[n] | mu, sigma)));
    
  }
  
}

model {
  
  y ~ normal(mu, sigma); // (0)
  
  for (n in 1:N) {
    
    // Jacobian adjustment
    target += log(sigma); // (1)
    
    target += log_diff_exp(normal_lcdf(y_star[n] + 1 | mu, sigma), 
                           normal_lcdf(y_star[n] | mu, sigma)); // (2)

    target += -normal_lpdf(y[n] | mu, sigma); // (3)
    
  }
  
}

```

This model doesn’t work well. Line `(2)` is the same as in the integrated case. Line `(3)` cancels out Line `(0)`. So the problem is Line `(1)` (the model works if I remove it), but I am not sure why.

Does my logic make sense? Where have I gone wrong?

---

<div class="post-metadata">

**Author:** ![aakhmetz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aakhmetz/32/18672_2.png) [@aakhmetz](https://discourse.mc-stan.org/u/aakhmetz)\
**Post date:** [July 2, 2025, 4:22am UTC](https://discourse.mc-stan.org/t/interval-censored-normal-model-integration-vs-estimation/39817/2 "2025-07-02T04:22:14Z")

</div>

I think you just introduce a new parameter, `y_raw` between 0 and 1, `y = y_raw + y_star[n]` with prior `y_raw ~ beta(1, 1)`. Then you fit it to your data such as (0). Maybe if this simpler model works you can move to more complicated variations

---

<div class="post-metadata">

**Author:** ![mdonoghoe](https://avatars.discourse-cdn.com/v4/letter/m/8dc957/32.png) [@mdonoghoe](https://discourse.mc-stan.org/u/mdonoghoe)\
**Post date:** [July 9, 2025, 5:45am UTC](https://discourse.mc-stan.org/t/interval-censored-normal-model-integration-vs-estimation/39817/3 "2025-07-09T05:45:59Z")

</div>

Thanks, I assume you meant that `y_raw` is a vector so `y[n] = y_raw[n] + y_star[n]`? That worked and is a nice approach that I hadn’t thought of, and although it doesn’t fit exactly into the more complex model that I’m working towards, it has given me some ideas.

In the meantime, I found the error I made above. My derivation of the Jacobian was correct, but I translated the -\log \phi((y\_i - \mu) / \sigma) as `-normal_lpdf(y[n] | mu, sigma)` in line `(3)` when it should have been `-std_normal_lpdf((y[n] - mu) / sigma)` , or equivalently `-normal_lpdf(y[n] | mu, sigma) - log(sigma)`, which would have cancelled out line `(1)`.