Hello there,

I need your help guys :S. I am working on a pretty simple linear model on real and imaginary components for recorded sound pressure \tilde{p} = \hat{p}+n where n is the noise. The model is a truncated sum of the spherical harmonics basis

where K and C are given constants. I want to find both A_{m,r} and A_{m,i} til some truncated M. I am modeling n\sim \mathcal{N}(0,1) and the A's are independent, modeled as A\sim \mathcal{N}(0,4). Here the Stan code

```
data {
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_meas] pt_real; // Measured Pressure at receiver. Real part
vector[N_meas] pt_imag; // Measured Pressure at receiver. Imaginary part
vector[N_meas] legendre[M]; // Legendre polynomials
vector[N_meas] bessel[M]; // Spherical bessel functions
vector[N_meas] neumann[M]; // Spherical Neumann functions
real<lower=0> e_real; // Prior variance from noise
real<lower=0> e_imag; // Prior variance from noise
real<lower=0> a_real; // variance of normal distributed A's
real<lower=0> a_imag; // variance of normal distributed A's
}
parameters {
vector[M] Ar; // real part of the source strength
vector[M] Ai; // imaginary part of the source strength
}
transformed parameters{
vector[N_meas] mu_real; // mean of the real part of the pressure
vector[N_meas] mu_imag; // mean of the imaginary part of the pressure
// Loop over modes
mu_real = 0 * d;
mu_imag = 0 * d;
for (m in 1:M){
mu_real += d .* (Ar[m] * bessel[m] + Ai[m] * neumann[m]) .* legendre[m];
mu_imag += d .* (Ai[m] * bessel[m] - Ar[m] * neumann[m]) .* legendre[m];
}
}
model {
Ar ~ normal(0, a_real);
Ai ~ normal(0, a_imag);
pt_real ~ normal(mu_real, e_real);
pt_imag ~ normal(mu_imag, e_imag);
}
```

This problem is solved for each frequency independently and the data has very little error (recorded in anechoic conditions). For every frequency I am experiencing the same problem when running inferences. Up to certain M the solution converges nicely but getting higher in M it makes it divergent (the M where it gets messed up depends on the frequency). It seems that it has problems having many parameters which should shrink towards 0. The amount of data is around 100-200 observations and the number modes M never goes higher than 10 so the problem is quite overdetermined.

I have solved the posterior analytically and everything works flawlessly, but it would be great if it works in Stan so I can change prior, hyperparameterize… The Stan output is the same as the analytic posterior until M is high enough so it diverges.

I upload the model and the data for a single frequency (100 Hz). For this case it converges up to M=5.

Thanks!

sh_simple.stan (1.8 KB)

data.csv (8.7 KB)