Slow, underlying trend with a random walk or Gaussian process

Question

I'm trying to fit a PyMC3 model to some data regarding sales over time. Here's a brief description :

• N salespeople each sell some number of widgets per week
• We assume each salesperson sells widgets at a different mean rate per week, and call this beta_i for salesperson i
• Our observed data is assumed to be ~Poisson(beta_i).

Weekly average sales data is plotted here in a histogram, with a log-normal fit on top, to give you an idea of the distribution of weekly widget sales by salesperson.

In this first scenario, I get what I think are a reasonable set of betas, although they don't look particularly log-normal :

Because we are hoping to infer something about an underlying trend shared by all salespeople (something analogous to "the economy"), we tried adding something. Our first attempt had "the economy" be just a linear function of time, starting at an intercept value of 1 and having derivative gamma > 0 ( gamma was half-normal with sd=0.5 ). We then had our data ~Poisson(beta_i * (1 + gamma)). In this scenario, betas didn't shift much, and we did infer something about "the economy", though it was a pretty weak effect.

I'm hoping to replace this with a random walk or a Gaussian process, to allow "the economy" to vary somewhat smoothly in time but to have an arbitrary shape. Ideally it would start at a value of 0, and then go wherever it needs to go to capture the underlying trend shared by all salespeople, with the data once again ~Poisson(beta_i * (1 + gamma)). Here's our model.

with pm.Model() as model:

    # Salesperson base rate of selling widgets
    beta_ = pm.Lognormal("beta", mu=mu_hat, sd=sd_hat, shape=(1, n_salespeople))
    # mu_hat and sd_hat were estimated by fitting a log-normal to weekly sales data

    # Economy
    gamma_ = pm.GaussianRandomWalk("gamma", mu=0, sd=1e-6, shape=(n_weeks, 1))


    # Effects
    base_rate = beta_
    economy = 1 + gamma_

    # Observed
    lambda_ = base_rate * economy
    y = pm.Poisson("y", mu=lambda_, observed=observed_sales + 1e-7)

where observed_sales is an integer array of the number of sales made, of shape (n_weeks, n_salespeople).

First, I'm not sure I'm specifying this model correctly. Without the "economy", I infer a reasonable set of betas ( although it doesn't look log-normal, as in the second screenshot ). Then, the random walk we get back is not at all smooth, no matter how small the sd gets; more often than not, for reasons I'm unsure about, I get a message about "Mass matrix contains zeros on the diagonal.". Finally, even at the beginning, I was getting infinite probabilities if I didn't add a small factor to the observed data... Why is that ?

So, TL; DR : I'm fairly new to probabilistic programming, and I'm fairly sure something is going wrong but I'm not sure what. Any input much, much appreciated !