[ENH] Vectorization support for zero-dimensional distributions

[tingiskhan]

Is your feature request related to a problem? Please describe.

Currently, zero-dimensional distributions don't support vectorization for methods s.a. log_pmf, pdf etc. Supporting this would make sense since this is supported by other similar libraries (like pytorch/scipy).

As an example, we'd go from this

from skpro.distributions import Normal
import numpy as np

dist = Normal(0, 1)
x = np.linspace(-2, 2)

pdf = np.array([dist.pdf(xt) for xt in x])

to

from skpro.distributions import Normal
import numpy as np

dist = Normal(0, 1)
x = np.linspace(-2, 2)

pdf = dist.pdf(x)

Describe the solution you'd like

The reason why this doesn't work is because of the below function

    def _coerce_to_self_index_np(self, x, flatten=False):
        """Coerce input to type similar to self.

        Coerces x to a np.ndarray with same shape as self.
        Broadcasts x to self.shape, if necessary, via np.broadcast_to.

        Parameters
        ----------
        x : array-like, np.ndarray coercible
            input to be coerced to self
        flatten : bool, optional, default=True
            if True, flattens x before broadcasting
            if False, broadcasts x as is

        Returns
        -------
        np.ndarray of same shape as self
            coerced input, with same index and columns as self, or float
            Return is ``x`` broadcasted to ``self.shape``, via ``np.broadcast_to``.
            If ``flatten`` is True, flattens ``x`` before broadcasting,
        """
        x = np.array(x)
        if flatten:
            x = x.reshape(1, -1)
        df_shape = self.shape
        x = np.broadcast_to(x, df_shape)
        return x

When the distribution is zero-dimensional, df_shape is an empty tuple and the broadcasting thus fails. I think it would make sense to skip the broadcasting if self is empty.

This also relates to #622.

[Aymuos22]

I'll soon raise a PR for the below:

I changed the internal coercion helpers so scalar distributions don’t try to broadcast every input to shape (), which means Normal(0, 1).pdf(np.linspace(...)) works without a Python loop. I also tightened the boilerplate so it only collapses true scalar outputs. Added regression tests for scalar Normal and Poisson that compare the vectorised pdf/log-pdf/cdf (or pmf/log-pmf) results against an element-wise baseline to ensure we don’t regress

[Aymuos22]

Hi @tingiskhan @fkiraly

Kindly check:
#676