Description of your problem or feature request

It is my understanding that current exact posteriors, e.g. gamma_poisson_conjugateo, can only condition on a single observation $Y \sim \text{Poisson}$ rather than a set of i.i.d. observations $Y_1, \dots, Y_n$. How can current conjugate relations be extended for multiple observations?

Can we allow arguments such as realized (akin to joint_logprob in AePPL) or realized_rvs_to_values in the AeMCMC's nuts' construct_sampler in AeMCMC's general construct_sampler?

srng = at.random.RandomStream(0)

lam_rv = srng.gamma(1., 1., name="lam")
Y_rv = srng.poisson(lam=lam_rv, size=3, name="Y") # something like this?

y_vv = Y_rv.clone()

sampler, initial_values = aemcmc.construct_sampler({Y_rv: y_vv}, srng)

p_posterior_step = sampler.sample_steps[lam_rv]

sample_fn = aesara.function([y_vv], p_posterior_step)

Currently, sample_fn(np.array([2, 3, 4]) yields a 3-dimensional array.

cc @brandonwillard @rlouf