Use scan to broadcast models fitting along the data axis

[LegrandNico]

This for loop is slowing model fitting when many datasets/models are provided.

https://github.com/ilabcode/ghgf/blob/8fec5d71e1677e85103ae9dd79dda26056b5b06f/ghgf/distribution.py#L167

We should use scan or for loop instead, but this might require careful refactoring.

It is not clear however if we would get significant performance improvement from this refactoring, if so it might only appear when n is large. Also, it might require uniformizing the shape of the input arrays, which will also come with a cost.