How to model sinusoids in white gaussian noise

Oh, and if you have frequencies with very different orders of magnitude, then maybe it is better to use VanMiseslength with some arbitraty length distribution that the sampler likes (say pm.Normal.dist(mu=4, sd=1) or so), then just ignore the second output and introduce a separate variable for the actual frequency:

with pm.Model() as model:
    length = pm.Bound(pm.Normal, lower=0).dist(mu=4, sd=1)  # The sampler seems to like this
    vals = VonMisesLength('a', mu=np.pi/4, kappa=1, length_dist=length, shape=(3, 2))
    # _ = vals[..., 0] # we don't need this
    phase = vals[..., 1]  # is in (-pi, pi)
    frequency = pm.HalfNormal('freq', ...)