def blend_means_sigmas(means, sigmas, weights):
    """Calculate the blended means and sigmas, the normalization parameters
    needed to recompose the cascade. This procedure uses the weights of the
    blending of the normalized cascades and follows eq. 32 and 33 in BPS2004.

    Parameters
    ----------
    means : array-like
      Array of shape [number_components, scale_level, ...]
      with the mean for each component (NWP, nowcasts, noise).
    sigmas : array-like
      Array of shape [number_components, scale_level, ...]
      with the standard deviation for each component.
    weights : array-like
      An array of shape [number_components + 1, scale_level, ...]
      containing the weights to be used in this routine
      for each component plus noise, scale level, and optionally [y, x]
      dimensions, obtained by calling either
      :py:func:`pysteps.blending.steps.calculate_weights_bps` or
      :py:func:`pysteps.blending.steps.calculate_weights_spn`.

    Returns
    -------
    combined_means : array-like
      An array of shape [scale_level, ...]
      containing per scale level (cascade) the weighted combination of
      means from multiple components (NWP, nowcasts and noise).
    combined_sigmas : array-like
      An array of shape [scale_level, ...]
      similar to combined_means, but containing the standard deviations.

    """
    # Check if the dimensions are the same
    diff_dims = weights.ndim - means.ndim
    if diff_dims:
        for i in range(diff_dims):
            means = np.expand_dims(means, axis=means.ndim)
    diff_dims = weights.ndim - sigmas.ndim
    if diff_dims:
        for i in range(diff_dims):
            sigmas = np.expand_dims(sigmas, axis=sigmas.ndim)
    # Weight should have one component more (the noise component) than the
    # means and sigmas. Check this
    if (
        weights.shape[0] - means.shape[0] != 1
        or weights.shape[0] - sigmas.shape[0] != 1
    ):
        raise ValueError(
            "The weights array does not have one (noise) component more than mu and sigma"
        )
    else:
        # Throw away the last component, which is the noise component
        weights = weights[:-1]

    # Combine (blend) the means and sigmas
    combined_means = np.zeros(weights.shape[1])
    combined_sigmas = np.zeros(weights.shape[1])
    total_weight = np.sum((weights), axis=0)
    for i in range(weights.shape[0]):
        combined_means += (weights[i] / total_weight) * means[i]
        combined_sigmas += (weights[i] / total_weight) * sigmas[i]

    return combined_means, combined_sigmas
Calculate the blended means and sigmas, the normalization parameters
needed to recompose the cascade. This procedure uses the weights of the
blending of the normalized cascades and follows eq. 32 and 33 in BPS2004.

Parameters
----------
means : array-like
Array of shape [number_components, scale_level, ...]
with the mean for each component (NWP, nowcasts, noise).
sigmas : array-like
Array of shape [number_components, scale_level, ...]
with the standard deviation for each component.
weights : array-like
An array of shape [number_components + 1, scale_level, ...]
containing the weights to be used in this routine
for each component plus noise, scale level, and optionally [y, x]
dimensions, obtained by calling either
:py:func:`pysteps.blending.steps.calculate_weights_bps` or
:py:func:`pysteps.blending.steps.calculate_weights_spn`.

Returns
-------
combined_means : array-like
An array of shape [scale_level, ...]
containing per scale level (cascade) the weighted combination of
means from multiple components (NWP, nowcasts and noise).
combined_sigmas : array-like
An array of shape [scale_level, ...]
similar to combined_means, but containing the standard deviations.
pysteps/blending/steps.py