I have some binned data (no access to underlying info) and prior knowledge that the value in each bin smoothly varies in space. So I am modeling using a Gaussian Process prior, which according to Rasmussen has the following marginal likelihood:

$$ p(\mathbf{y}) = (2\pi)^{-N/2}{|K|}^{-1/2}\exp \left(-\frac{1}{2}(\mathbf{y}-\mathbf{m})^\intercal K^{-1}(\mathbf{y}-\mathbf{m})\right) $$ Where $K$ is the covariance, $y$ are the observations, and $m$ is the mean function. How should we adjust this when instead of point observations $\mathbf{y}$ we have average values for finite volumes (in 3D)?

What I have so far: Perhaps the answer is to replace $\mathbf{y}-\mathbf{m}$ with the difference vector with entries: $$ \begin{align} d_i &= \int_{V_i} (y_i-m(\mathbf{r}))d\mathbf{r} \\ &\approx |V_i|(y_i-m_i) \end{align} $$

Where the last approximation has volume $|V_i|$ and $m_i$ is the mean function evaluated at the bin center.

A related question: how would we sample densities from this distribution?

  • $\begingroup$ What is the dimension of the sample space? I.e. are these bivariate normals or higher (and how much higher)? $\endgroup$
    – Cliff AB
    Commented Mar 25, 2016 at 0:09
  • $\begingroup$ @CliffAB The samples are embedded in 3D space. $\endgroup$
    – cgreen
    Commented Mar 25, 2016 at 0:12
  • $\begingroup$ Just curious: is this for an engineering-based application by any slim chance? I talked with some engineers who stated what sounded like a very similar problem, but haven't yet followed up with them on that. But maybe that's just a coincidence. $\endgroup$
    – Cliff AB
    Commented Mar 25, 2016 at 0:56
  • $\begingroup$ @CliffAB thanks, this is a biology application actually. Interpreting microscopy data. $\endgroup$
    – cgreen
    Commented Mar 25, 2016 at 3:42

1 Answer 1


Frequent readers of Cross-Validated must be beginning to think that I'm some sort of broken record at this point, but this can be considered interval censored data. In other words, because of the binning, we know that each value (in each dimension) lies in a given interval, but not the exact value. See this post for a very similar problem (the same problem, but only in 1 dimension).

This is a little trickier, because of the higher dimensions. In this case, let us define $L_i^j$ to be the lower limit of the $j^{th}$ dimension for the $i^{th}$ subject and $U_i^j$ be the corresponding upper limit of the interval. As an example, suppose our binning is from simple rounding and the unobserved tuples are (1.2, 0.2, -0.5). After binning, we only know that is the original values are in the following intervals: [1.0, 2.0), [0.0, 1.0) and [-1.0, 0.0). If this were for the first subject, we would then say $L_1^1 = 1.0$, $U_1^1 = 2.0$, $L_1^2 = 0.0$, $U_1^2 = 1.0$, etc.

Now, if $f(x_1, x_2, x_3 | \Theta)$ is the multivariate density function for the joint distribution of the exact process (i.e. the true value before binning), then we can write the likelihood function for the binned data as

$L(\Theta | \text{data}) = \displaystyle \prod_{i = 1}^n \int_{L_i^1}^{U_i^1} \int_{L_i^2}^{U_i^2} \int_{L_i^3}^{U_i^3} f(x_1, x_2, x_3 | \Theta) dx_3 dx_2 dx_1$

(note: this assumes that each dimension is binned. If some of the dimensions are binned and some are not, it gets a little nastier.)

From here, standard MLE or Bayesian methods can be used for inference. If you want to fit the Gaussian model, you just use the multivariate Gaussian function.

However, easier said than done! I'm not aware of any statistical packages that are built specifically for fitting this model. If you don't mind rolling up your sleeves and handwriting the model, you should be able to fit these models in Stan, JAGs or nimble. Note if you go that route, you can use their syntax for interval censoring and do not need to write out the integrals.

Alternatively, you could use a multiple imputation based methods. I don't believe there's any currently available software packages for that particular model.

The good news is that once you define the underlying distribution (i.e. the distribution of the un-binned points), sampling of the binned points is done very easily. Just first sample the un-binned points and then place them in the correct bin, as I did with my simple rounding binning example.


So if we have a current value of $\mu$ and $\Sigma$, how do we calculate the integral for the likelihood function? Well, let's assume we the multivariate normal CDF function available (for example, this is available in R's mvtnorm package through the function pmvnorm). I will denote this function as $\Phi(x_1, x_2, x_3, ...| \mu, \Sigma)$, (and assume $\mu$, $\sigma$ are the correct corresponding dimensions as well).

Computing the integral for 1D is straight forward: this is just

$ \Phi(R^1_i | \mu, \Sigma) - \Phi(L^1_i| \mu, \Sigma)$

For 2D, it gets a little hairier:

$ \Phi(R^1_i, R^2_i | \mu, \Sigma) - \Phi(R^1_i, L^2_i | \mu, \Sigma) - \Phi(L^1_i, R^2_i | \mu, \Sigma) + \Phi(L^1_i, L^2_i | \mu, \Sigma)$

For 3D...well, thank goodness your not going past 3D. Here's 3D:

$ \Phi(R^1_i, R^2_i, R^3_i | \mu, \Sigma) - \Phi(L^1_i, R^2_i, R^3_i | \mu, \Sigma) -\Phi(R^1_i, L^2_i, R^3_i | \mu, \Sigma) -\Phi(R^1_i, R^2_i, L^3_i | \mu, \Sigma) +\Phi(L^1_i, L^2_i, R^3_i | \mu, \Sigma) +\Phi(L^1_i, R^2_i, L^3_i | \mu, \Sigma) + \Phi(R^1_i, L^2_i, L^3_i | \mu, \Sigma) - \Phi(L^1_i, L^2_i, L^3_i | \mu, \Sigma)$

  • $\begingroup$ Also, note that the stated likelihood function is equivalent to finding the total probability mass of the cube dictated by the 3D bin. So I believe this is approximately what you were stating about finding the volume (over the density function) of the bin containing each point. $\endgroup$
    – Cliff AB
    Commented Mar 25, 2016 at 1:12
  • $\begingroup$ thanks so much for your response! Makes a lot of sense. One issue is that my voxels aren't independent, so the total likelihood is not just the product over each voxel likelihood, but rather has covariance (seen in my Gaussian process likelihood). So, basically, I don't know where to put the integral. $\endgroup$
    – cgreen
    Commented Mar 25, 2016 at 3:45
  • $\begingroup$ @cgreen: I wrote out explicitly how one would calculate the integral (well, how you would calculate it with available software, anyways). $\endgroup$
    – Cliff AB
    Commented Mar 25, 2016 at 4:19
  • $\begingroup$ Thanks, however, I think what you have written is the likelihood of a gaussian density on a histogram. What I need is the likelihood of a gaussian process, that is, a collection of random variables with a joint multivariate distribution. That is, each point is linked to each other point via the covariance $\Sigma$, so the likelihood is not $\prod_i N(y_i|m,\sigma)$ but rather the multivariate $N(\mathbf{y}|\mathbf{m},\Sigma)$. The problem I have is that the variables are spread out over a bin. $\endgroup$
    – cgreen
    Commented Mar 25, 2016 at 5:19
  • $\begingroup$ @cgreen: Oh okay, so the last 3D time point is correlated with the next 3D time point, like a state-space model? Oh man, that likelihood is going to get uncomfortable quickly. You should still be able to do this in the BUGs framework (i.e. Stan, JAGs, etc), in which case writing out the likelihood won't be so bad and you can use the generic samplers but I would be willing to bet that the mixing of your samples will be terrible. How many samples are we talking about? $\endgroup$
    – Cliff AB
    Commented Mar 25, 2016 at 5:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.