Thanks a lot for promptly replying.
Sorry, it’s unclear, but that’s the point I couldn’t connect the model to implementation… I’ll try to elaborate on that. Actually, I intend to do something like this
X \thicksim \textsf{Uniform}(0, L)
Y \thicksim \textsf{Uniform}(0, W)
D_i \thicksim \sqrt{\left( X - x_i \right)^2 + \left( Y - y_i \right)^2}
\mu_i \thicksim \alpha + \beta \log(D_i)
\alpha \thicksim \textsf{Normal}(0, 100)
\beta \thicksim \textsf{Normal}(0, 100)
\sigma^2 \thicksim \textsf{HalfNormal}(10)
\rho_i \thicksim \textsf{Normal}(\mu_i, \sigma^2)
So that, inside the sampling algorithm over the updating iterations (I guess)
\displaystyle X_{(k + 1)} \thicksim \Pr(X)\prod_{i=1}^N \Pr \left( \rho_i \vert Y_{(k)}, \sigma^2_{(k)}, \alpha_{(k)}, \beta_{(k)} \right)
\displaystyle Y_{(k + 1)} \thicksim \Pr(Y)\prod_{i=1}^N \Pr \left( \rho_i \vert X_{(k + 1)}, \sigma^2_{(k)}, \alpha_{(k)}, \beta_{(k)} \right)
Then, the corresponding conditional distributions would be written as
\displaystyle X_{(k + 1)} \thicksim \exp\{-\frac{\sum_{i=1}^N \left[ \rho_i - \alpha_{(k)} - \beta_{(k)} \log(D_i)\right]^2}{2 \sigma^2_{(k)}}\}, \text{with } D_i \thicksim \sqrt{\left( X - x_i \right)^2 + \left( Y_{(k)} - y_i \right)^2}
\displaystyle Y_{(k + 1)} \thicksim \exp\{-\frac{\sum_{i=1}^N \left[ \rho_i - \alpha_{(k)} - \beta_{(k)} \log(D_i)\right]^2}{2 \sigma^2_{(k)}}\}, \text{with } D_i \thicksim \sqrt{\left( X_{(k + 1)} - x_i \right)^2 + \left( Y - y_i \right)^2}
I tried to implement that as a multilevel model (as discussed in A Primer on Bayesian Methods for Multilevel Modeling), but it seems to be stuck to the priors because I always get (50, 50) as an estimate to (X, Y) when considering L = 100\! m and W = 100\! m.