Preserve model parameter axes in censored Gaussian gradients

[antonsoo]

CensoredGaussianLogLikelihood.evaluateS1 assumes the number of model parameters equals the number of outputs. For example, the two-output, four-parameter LotkaVolterraModel raises a broadcast error. With matching counts, parameters shared by outputs can instead receive incorrect gradients because the censored terms weight and sum the wrong axis.

Broadcast the censoring masks over model parameters, apply each output's noise scale, and sum over outputs. Keep the likelihood values and model evaluations unchanged.

Add 24 numerical regression cases covering one/two outputs, one/two/three model parameters, shared parameters, unequal noise, and no/lower/upper/both censoring. Seventeen cases fail on unchanged main; all pass after on Python3.10/NumPy1.26 and Python3.12/NumPy2.5.

I also checked 24 continued gradient-descent fits using the recorded hare/lynx series supplied with LotkaVolterraModel, constructed censoring variants, and ordinary Gaussian/diagonal controls. Twelve shape errors and three incorrect gradient trajectories become zero against an independent ODE solve and finite-difference reference; all nine previously working control trajectories stay identical.

Validation: full unit suite passes (766 tests, three optional skips), Flake8 passes, and the clean HTML documentation build succeeds without warnings. CHANGELOG.md is updated.

Related to #667.