Gaussian process likelihood with binned data

Question

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?

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.

EDIT

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)$