Solving independent problems (frequencies) at once

Modeling
1

Hello everyone,

I am working on soundfield reconstruction. I am using an orthogonal basis (spherical harmonics) to model some measurements from one loudspeakers to several microphones. The basis is the multiplication of a radial component (spherical bessel and neumann functions) and an angular component (legendre polynomials). I want to find the coefficients of the different modes in the basis (named in the code Ar and Ai for real and imaginary).

data {
  int<lower=0> N;                             // Number of frequencies.
  int<lower=0> N_meas;                        // Total Number of independent measurements.
  int<lower=0> M;                             // Number of spherical harmonics modes.                                   
  vector[N_meas] d;                           // Distance loudspeaker-microphone.
  vector[N] pt_real[N_meas];                  // Measured Pressure at receiver. Real part
  vector[N] pt_imag[N_meas];                  // Measured Pressure at receiver. Imaginary part
  real legendre[M, N_meas];                   // Legendre polynomials
  vector[N] bessel[M, N_meas];                // Spherical bessel functions
  vector[N] neumann[M, N_meas];               // Spherical Neumann functions
}
parameters {
  vector<lower=0>[N] epsilon_real;            // variance of the background noise. Real part.
  vector<lower=0>[N] epsilon_imag;            // variance of the background noise. Imaginary part.
  vector[N] Ar[M];                            // real part of the source strength
  vector[N] Ai[M];                            // imaginary part of the source strength  
}
transformed parameters{
  vector[N] mu_real[N_meas];                  // mean of the real part of the pressure
  vector[N] mu_imag[N_meas];                  // mean of the imaginary part of the pressure
  // Loop over frequencies
  for (n in 1:N){
    // Loop over measurements
    for (i in 1:N_meas){
      mu_real[i][n] = 0;
      mu_imag[i][n] = 0;
      for (m in 1:M){
        mu_real[i][n] += d[i] * (Ar[m][n] * bessel[m, i][n] + Ai[m][n] * neumann[m, i][n]) * legendre[m, i];
        mu_imag[i][n] += d[i] * (Ai[m][n] * bessel[m, i][n] - Ar[m][n] * neumann[m, i][n]) * legendre[m, i];
      }
    }
  }
}
model {
  epsilon_real ~ normal(0, 1);
  epsilon_imag ~ normal(0, 1);
  for (m in 1:M){
    Ar[m] ~ normal(0, 4);
    Ai[m] ~ normal(0, 4);
  }
  for (i in 1:N_meas){
    pt_real[i] ~ normal(mu_real[i], epsilon_real);
    pt_imag[i] ~ normal(mu_imag[i], epsilon_imag);
  }
}

The problem is iid and independent per frequency (dimension N in the code), so I would expect that solving for [100] Hz frequency or for [100, 200] Hz at the same time would give me the “same” result for Ar and Ai at 100 Hz. However, when I introduce lower frequencies (which are more difficult to fit) it messes up all the results, including other frequencies (It doesn’t converge anymore). Probably my code is wrong or my understanding on how Stan works is wrong, or both. Do you have any answer to this?

Thanks!