McCullagh_compute_cumulatives function

Computes the model-based cumulative probability matrices pij and qij