Because the above model runs surprisingly slow, I decided to simplify the model further to be
amplt = pm.Normal(‘amplt’, mu=1, sd=1)
freq = pm.Uniform(‘freq’, 0, 0.5)
mu_x = amplt * pm.math.sin(2*pi * t * freq)
X_obs = pm.Normal(‘X_obs’, mu=mu_x, sd=1, observed=x)
and it still took 18:40 to take 5000 samples with NUTS. The number of data points is 512. After about 2500 samples, it hit on the correct answer pair. I presume that if I ran the simulation longer, the posterior samples would likely be far more weighted to
I guess the underlying question is how can I speed up this inference to something faster than ~4.5 iterations per second. Lots of important problems aren’t a single 0-deg phase sinusoid but a complex frequency structure that may require a rather messy multi-level model. I expect there are methods one can leverage to improve performance but I am still somewhat new to this API.