Hi all.

I’m trying to fit multi-level wave regressions for my research. Originally I was trying to do this with brms in R but got some guidance that I might want to try a library that supports SMC, because the posterior is multimodal and SMC is better suited for this geometry. I moved to PyMC.

Currently I’m working on simulated data. It is better than what I was getting previously with brms’ NUTS implementation though still not really converging in the way I would hope to see. I was interested to get guidance from the forums here.

Here’s a sample of the kind of thing I’ve been trying to fit. I based this approach on the pymc multi-level model tutorial I found here:

```
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
# generate wave data
np.random.seed(1)
params_dict = {
"ID": list(range(100)),
"freq": np.random.normal(1/7, 0.05, 100),
}
df_params = pd.DataFrame(params_dict)
ids = np.repeat(df_params['ID'].values, 20)
timepoints = np.tile(np.arange(0,20), len(df_params))
df_new = pd.DataFrame({'ID': ids, 'time': timepoints})
# Step 2: Merge with original DataFrame
df_data = pd.merge(df_new, df_params, on='ID')
epsilon = np.random.normal(0, 0.1, len(df_data))
df_data['y'] = np.cos(2 * np.pi * df_data['freq'] * df_data['time']) + epsilon
# Step 3: Create model
coords = {"part": df_data["ID"].unique()}
part = df_data["ID"].values
y = df_data["y"].values
with pm.Model(coords=coords) as model:
x = pm.MutableData("x", df_data["time"], dims="obs_id")
part_idx = pm.MutableData("part_idx", part, dims="obs_id")
freq_fixed = pm.Normal("freq_fixed", mu=1/7, sigma=0.05)
freq_sd = pm.Normal("freq_sd", 0.05, 0.01)
freq = pm.Normal("freq", mu=freq_fixed, sigma=freq_sd, dims="part")
eps = pm.Normal("eps", 0.1, 0.01)
mu = np.cos(2*np.pi*freq[part_idx]*x)
y_like = pm.Normal("y_like", mu=mu, sigma=eps, observed=y, dims="obs_id")
# Step 4: Sample
with model:
trace = pm.sample_smc(2000,
random_seed=4,
progressbar=False,
cores=4)
```

So my data look something like this:

And the model graph looks like this:

model_graph.pdf (3.7 KB)

The fits for this model are better than what I was getting when I was using HMC, but still definitely not fully converged. Rhats are consistently above 1.01. Here’s a sample trace plot:

One thing that is kind of interesting is that the parameter values it seems to be finding are not necessarily the ones I used to generate the data. It seems to be getting pretty close to the frequency mean, but the standard deviation of the frequency is usually around 0.04 and the true value is 0.05. Similarly the epsilon value is often inferred around .2 when the true value is 0.1. This is in spite of the fact that I’m specifying pretty strong priors on the true values that I used to generate the data.

When I weaken the priors or try to add additional wave parameters to infer simultaneously performance starts to deteriorate. I’m wondering whether there is anything I can do to improve things here? Or perhaps other algorithms I should try? Appreciate any guidance! Thank you!