I have a hierarchical model that I’m not sure how to deal with and wanted to ask for some tips.

Below is a graphical model showing the gist. I would like to learn about alpha and beta. There is an individual per-object parameter y_i that is related to alpha and beta, and another parameter z_i that is unrelated to alpha and beta. By related, I mean the relationship between x_ij and y_i is dependent on alpha and beta.

The problem is that my model for how alpha, beta, and y_i generates x_ij is not analytic and expensive to compute. I only have it for certain values and need to interpolate between values (based on the assumption that it should vary continuously).

Once alpha and beta are fixed, I can approximate my model as e.g., spline and I think I can code up the things inside j plate in stan. Any thoughts would be appreciated.

Hi, this is IMHO very hard to understand without knowing what the model for how \alpha, \beta, and y_i generates x_{ij} looks like. Generally Stan requires a generative model (a model you could run to simulate new data from known parameter values) so if you can’t express this, you would likely need an analytic approximation or you are in trouble.
You could in principle let Stan fit the function f: \{\alpha,\beta,y_i\} \rightarrow \{x_{ij}\} as a trivariate spline/Gaussian process, but AFAIK 3D Gaussian processes are very hard to fit and splines might not be better + the problem may not really be identified.

Thanks for the reply. I realized that because f needs to be interpolated from three-dimensional grid (in which each point is costly to compute), I think I need to write this as an external c++ function and provide that to stan, which seems like a difficult task.

I didn’t know Stan can do multidimensional spline interpolation internally. Could you direct me to where I can read about that? (or maybe you meant do that with stan maths myself)

Stan AFAIK can’t do splines internally - you would have to implement that yourself, but it should be possible to do it in Stan, without resorting to C++. I honestly don’t even really know how multidimensional splines work, but a quick Google turned up some descriptions that look feasible at first glance. It’s likely not gonna be simple, but likely feasible and not extremely difficult.

I only have experience with univariate splines in Stan and those are +/- easy and code is available so you might be able to start off of that (https://github.com/milkha/Splines_in_Stan). There is also a Stan example with sparse splines representation (where most coefficients are 0) which might be helpful if you end up with a lot of splines (see Spline fitting demo inc comparison of sparse vs non-sparse)

Actually, I realized there are two cases which might have been confused:
a) You know that the values at the grid are exactly correct, and you need to just evaluate a spline, where the coefficients have been computed beforehand (possibly outside Stan)
b) The values at the grid are only approximate and you want to fit spline coefficients while fitting to the observed y_{ij}.

My above answer was mostly concerned with b), which is harder. If you can do the fit beforehand (outside Stan) and just need to evaluate the spline for given parameter values in the model block (case a), it should be easier - especially if a multidimensional spline fitting/interpolating package is available for your language. But in both cases you would need to implement the spline evaluation yourself in Stan.