Nothing urgent, but maybe someone has some comments on the set up.

In an initial experiment, proposed by @anon75146577 here (Algebraic sovler problems (Laplace approximation in Stan)), we had a Poission latent Gaussian model:

\phi \sim \mathrm{Normal}(0, \lambda) \\ \theta_i \sim \mathrm{Normal}(0, \Sigma_\phi) \\ y_{j \in g(i)} \sim \mathrm{Poisson}(e^{\theta_i})

where \Sigma_\phi was a diagonal covariance matrix, with ith entry \sigma_i = \phi^2. Now, as proposed by @avehtari, we want to generalize the problem a bit, and fit a model with a dense covariance matrix. So how do we generate this covariance?

My proposition is use as diagonal elements \sigma_{ii}^2 = \phi^2 + \epsilon, where \epsilon is a random perturbation, and for the off-diagonal elements use:

\sigma_{ij}^2 = \rho \times \sqrt{\sigma_i^2 \sigma_j^2}

where \rho is a correlation coefficient, which will vary during the experiment. Ok, should I add another perturbation to \sigma_{ij}? Is it ok for \Sigma_\phi to be a stochastic function of \phi?