The model is fitted to pairs of sites J and K independently (J, K = B1, B2, D1, D2), yielding a sample from the marginal distribution, e.g., for b J and b K (and similarly for the other parameters) for each of the sites respectively.