def create_P(P_δ, P_ζ, P_ι):
    """
    Combine `P_δ`, `P_ζ` and `P_ι` into a single matrix.

    Parameters
    ----------
    P_δ : ndarray(float, ndim=1)
        Probability distribution over the values of δ.

    P_ζ : ndarray(float, ndim=2)
        Markov transition matrix for ζ.

    P_ι : ndarray(float, ndim=1)
        Probability distribution over the values of ι.

    Returns
    ----------
    P : ndarray(float, ndim=3)
        Joint probability distribution over the values of δ, ζ and ι.
        Probabilities vary by δ on the first axis, by ζ on the second axis,
        and by ι on the third axis.

    """

    P = \
        P_δ[:, None, None, None] * P_ζ[None, :, :, None] * \
        P_ι[None, None, None, :]

    return P