For example, the parameter distribution recovered for the venue_split variable, which are 2 simple values, looks like this, so clearly the sampler has issues sampling this model, but I can’t tell where it’s misspecified:
Also, I’ve run a test model where I specify hard values for the alpha and beta parameters for each of the beta prior, which converges very well, so I’m thinking I’ve maybe misspecified the hyper-priors in the model above?
This is the test model, for comparison:
with pm.Model() as baseline_model:
# ---------------------------
# Reservation Rate
# ---------------------------
res_rate = pm.Beta("res_rate", alpha=res_rate_obs * 1000,
beta=(1-res_rate_obs)*1000,
shape=n_res_rate,
testval=res_rate_obs)
# ---------------------------
# Booking Rate
# ---------------------------
booking_rate = pm.Beta("booking_rate",
alpha=book_rate_obs * 1000,
beta=(1-book_rate_obs)*1000,
shape=n_event_weeks,
testval=book_rate_obs)
# ---------------------------
# Venue Split
# ---------------------------
venue_split = pm.Beta("venue_split",
alpha=venue_split_obs * 1000,
beta=(1-venue_split_obs)*1000,
shape=n_venues)
# ---------------------------
# Bookings
# ---------------------------
total_events = reservations_obs / res_rate[res_week_idx]
venue_bookings_rate = total_events * venue_split[venue_idx] * booking_rate[event_week_idx]
venue_bookings = pm.Poisson('venue_bookings', mu=venue_bookings_rate, observed=total_events_venue_obs)
This kind of convergence is what I’d like to see from the real model:

