def resample_distributions(
    first_array, second_array, probability_first_array, randgen=np.random
):
    """
    Merges two distributions (e.g., from the extrapolation nowcast and NWP in the blending module)
    to effectively combine two distributions for probability matching without losing extremes.
    Entries for which one array has a nan will not be included from the other array either.

    Parameters
    ----------
    first_array: array_like
        One of the two arrays from which the distribution should be sampled (e.g., the extrapolation
        cascade). It must be of the same shape as `second_array`. Input must not contain NaNs.
    second_array: array_like
        One of the two arrays from which the distribution should be sampled (e.g., the NWP (model)
        cascade). It must be of the same shape as `first_array`. Input must not contain NaNs.
    probability_first_array: float
        The weight that `first_array` should get (a value between 0 and 1). This determines the
        likelihood of selecting elements from `first_array` over `second_array`.
    randgen: numpy.random or numpy.RandomState
        The random number generator to be used for the binomial distribution. You can pass a seeded
        random state here for reproducibility. Default is numpy.random.

    Returns
    -------
    csort: array_like
        The combined output distribution. This is an array of the same shape as the input arrays,
        where each element is chosen from either `first_array` or `second_array` based on the specified
        probability, and then sorted in descending order.

    Raises
    ------
    ValueError
        If `first_array` and `second_array` do not have the same shape.
    """

    # Valide inputs
    if first_array.shape != second_array.shape:
        raise ValueError("first_array and second_array must have the same shape")
    probability_first_array = np.clip(probability_first_array, 0.0, 1.0)

    # Propagate the NaN values of the arrays to each other if there are any; convert to float to make sure this works.
    nanmask = np.isnan(first_array) | np.isnan(second_array)
    if np.any(nanmask):
        first_array = first_array.astype(float)
        first_array[nanmask] = np.nan
        second_array = second_array.astype(float)
        second_array[nanmask] = np.nan

    # Flatten and sort the arrays
    asort = np.sort(first_array, axis=None)[::-1]
    bsort = np.sort(second_array, axis=None)[::-1]
    n = asort.shape[0]

    # Resample the distributions
    idxsamples = randgen.binomial(1, probability_first_array, n).astype(bool)
    csort = np.where(idxsamples, asort, bsort)
    csort = np.sort(csort)[::-1]

    # Return the resampled array in descending order (starting with the nan values)
    return csort
Merges two distributions (e.g., from the extrapolation nowcast and NWP in the blending module)
to effectively combine two distributions for probability matching without losing extremes.
Entries for which one array has a nan will not be included from the other array either.

Parameters
----------
first_array: array_like
One of the two arrays from which the distribution should be sampled (e.g., the extrapolation
cascade). It must be of the same shape as `second_array`. Input must not contain NaNs.
second_array: array_like
One of the two arrays from which the distribution should be sampled (e.g., the NWP (model)
cascade). It must be of the same shape as `first_array`. Input must not contain NaNs.
probability_first_array: float
The weight that `first_array` should get (a value between 0 and 1). This determines the
likelihood of selecting elements from `first_array` over `second_array`.
randgen: numpy.random or numpy.RandomState
The random number generator to be used for the binomial distribution. You can pass a seeded
random state here for reproducibility. Default is numpy.random.

Returns
-------
csort: array_like
The combined output distribution. This is an array of the same shape as the input arrays,
where each element is chosen from either `first_array` or `second_array` based on the specified
probability, and then sorted in descending order.

Raises
------
ValueError
If `first_array` and `second_array` do not have the same shape.
pysteps/postprocessing/probmatching.py