Hello all!

Brand new to stan and HMC and having fun.

I have a spatial model of the form Y|X \sim N(|X|, \Sigma_2) and X \sim N(0, \Sigma_1), where X and Y are of length n and \Sigma_1, \Sigma_2 are fixed n \times n spatial covariance matrices. Upon a single observations of Y, the goal is to sample X|Y.

Here’s my code implementing the sampler:

```
data {
int n;
vector[n] Y;
matrix[n,n] Sigma1;
matrix[n,n] Sigma2;
vector[n] zeros;
}
parameters {
vector[n] X;
}
transformed parameters{
vector[n] absX = fabs(X);
}
model {
target += multi_normal_lpdf(X | zeros, Sigma1);
target += multi_normal_lpdf(Y | absX, Sigma2);
}
```

Evaluating the target distribution, a product of two multivariate-normals, requires solving linear systems with covariace matrices \Sigma_1 and \Sigma_2, which gets prohibitive for large n. In spatial statistics, we like to assume some computationally amenable form for the covariance matrices that allow for a fancy algorithm to efficiently solve the linear systems.

Is there a good way to incorporate a fancy solver into my stan implementation? Basically, I’d like the program to evaluate ‘multi_normal_pdf()’ in a special, efficient way. For example, could I write a seperate C++ function that stan would call to evaluate ‘my_special_multi_normal_pdf()’?

Thanks,

Paul