Apologies for the late reply. By using pm.sample_posterior_predictive I am only able to sample categorical outputs, over which we cannot compute values such as variance.
I did, however, manage to sample the probabilities from the model output using
sample_proba = approx.sample_node(model.out.distribution.p,
size=n,
more_replacements={input_var: x})