I’m trying to fit an image of (matrix of NxN) dimension, but I don’t figure out how to do it. I suppose that my model has to be a y ~multi_normal (mu, Sigma) with mu a vector of values mu_x and mu_y and Sigma a cov matrix.

I’m writing a simple Stan code, but i don’t know how to introduce the image to the multi_normal. An example of data could be:

[,1] [,2] [,3] [,4]

[1,] 1 4 7 6

[2,] 2 11 9 4

[3,] 3 6 9 4

[4,] 1 2 3 1

And I want to fit a 2d gaussian over this matrix. I’m lost… I would appreciate any clue.

data {
int<lower=0> n; #dimensions of gaussian
int<lower=0> N; #dimensions of the image
matrix[N,N] x; #image
}
parameters {
vector[2] mu;
matrix <lower=0> cov;
}
model {
y ~ multi_normal( mu, cov);
mu ~ normal(0., 100.); #prior
cov ~ cauchy(0., 100.);
}

As I’ll describe in a second this kind of model is really quite complicated and I recommend starting with something a bit simpler to help ease your transition into Stan. For example the cauchy density function is defined over one-dimensional variables but your Stan program is trying to evaluate it on the multi-dimensional variable cov. Similarly matrix<lower=0> is not a well-defined type in Stan.

I can come up with (at least) two ways of interpreting “fitting a 2d gaussian over a matrix”.

One is that you have a two-dimensional grid of pixels – a total of N^{2} pixels – and you want to model the value of each pixel with a joint distribution over all of those pixels. If those values can take on both negative and positive values then they can be modeled with a multivariate normal density function. This would require a mean vector with N^{2} elements, one for each pixel, and a covariance matrix with N^{2} \cdot N^{2} = N^{4} elements for each pair of pixels. Smoothness in the pixel values would correspond to strong correlations in neighboring pixels which would correspond to particular structure in the covariance matrix. Building covariance matrices (or precision matrices) with this kind of structure is a key focus of spatial statistical methods.

The other is that you want the pixels values to capture the image of a two-dimensional normal density function. Firstly what two-dimensional normal density function is being imaged? Do you want to fix a particular configuration (mu and cov in the data block) or is that configuration uncertain (mu and cov in the parameters block)? Once you’ve determined what is being imaged you need a model the location and extent of each pixel relative to the image and then how that image manifests in pixel values. For example you could integrate the normal density function across the spatial extent of each pixel (note that even in the two-dimensional scale these integrals cannot be done analytical) and then have the pixel value vary around that integral,

x_{n1, n2} \sim normal( I_{n1, n2}, sigma ).

Alternatively you could model the pixel values as concentrating around the value of the normal density function at the center of the pixel,

x_{n1, n2} \sim normal( y_{n1, n2}, sigma ).

In either case if you want to learn the latent image from observed pixel values then that requires another layer of complexity.

Anyways I know that this may not be the most satisfying answer but I hope it demonstrates that models are often much more complicated that we naively expect, but with careful though and consideration we can work out the relevant details. Moreover with experience with Stan we can implement those details in a Stan program.