This looks good to me. The sampling is well behaved and it reproduces your simulation values. It seems to me the questions you’re asking are more philosophical—“is it appropriate to call it ‘Poisson distributed’?” I suppose that depends on what you’re implying is your independent variable. Personally, I think of your data as Lorentz distributed with a Poisson variance. So, for a single x-value on your plot (or, practically speaking, for a narrow bin), you could look at the count distribution of the data in that bin and it’s reasonable to assume it’s Poisson distributed. But the value of \lambda for that Poisson will depend on your value of x (of course that corresponds to each value of \mu in your obs variable). Basically, in your code you’re fitting len(x) number of Poisson distributions with your obs variable.
Alternatively, you could convert your data into a 1D list and use the pm.Cauchy() distribution as your log-likelihood. For example say you had the following data points y=[(-10.5, 3), (-1, 11), (5.5, 4)] which would be re-expressed as y=[-10.5, -10.5, -10.5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 5.5, 5.5, 5.5, 5.5]. This 1D representation could then be set as the observed data for an explicitly Lorentz (a.k.a. “Cauchy”) likelihood.
i.e. obs = pm.Cauchy('obs', alpha=a_prior, beta=b_prior, observed=y)
I’m not sure if this is answering your question or if it’s even helpful, but perhaps it gives you another way of conceptualizing the problem. All of that said, I like the way you’re doing it.