def Sens_m_sample(poly, dist, samples, rule="R"):
    """
    First order sensitivity indices estimated using Saltelli's method.

    Args:
        poly (chaospy.Poly):
            If provided, evaluated samples through polynomials before returned.
        dist (chaopy.Dist):
            distribution to sample from.
        samples (int):
            The number of samples to draw for each matrix.
        rule (str):
            Scheme for generating random samples.

    Return:
        (numpy.ndarray):
            array with `shape == (len(dist), len(poly))` where `sens[dim][pol]`
            is the first sensitivity index for distribution dimensions `dim` and
            polynomial index `pol`.

    Examples:
        >>> dist = chaospy.Iid(chaospy.Uniform(), 2)
        >>> poly = chaospy.basis(2, 2, dim=2)
        >>> print(poly)
        [q0^2, q0q1, q1^2]
        >>> print(numpy.around(Sens_m_sample(poly, dist, 10000, rule="M"), 4))
        [[0.008  0.0026 0.    ]
         [0.     0.6464 2.1321]]
    """
    dim = len(dist)

    generator = Saltelli(dist, samples, poly, rule=rule)

    zeros = [0]*dim
    ones = [1]*dim
    index = [0]*dim

    variance = numpy.var(generator[zeros], -1)

    matrix_0 = generator[zeros]
    matrix_1 = generator[ones]
    mean = .5*(numpy.mean(matrix_1) + numpy.mean(matrix_0))
    matrix_0 -= mean
    matrix_1 -= mean

    out = [
        numpy.mean(matrix_1*((generator[index]-mean)-matrix_0), -1) /
        numpy.where(variance, variance, 1)
        for index in numpy.eye(dim, dtype=bool)
    ]

    return numpy.array(out)
First order sensitivity indices estimated using Saltelli's method.

    Args:
        poly (chaospy.Poly):
            If provided, evaluated samples through polynomials before returned.
        dist (chaopy.Dist):
            distribution to sample from.
        samples (int):
            The number of samples to draw for each matrix.
        rule (str):
            Scheme for generating random samples.

    Return:
        (numpy.ndarray):
            array with `shape == (len(dist), len(poly))` where `sens[dim][pol]`
            is the first sensitivity index for distribution dimensions `dim` and
            polynomial index `pol`.

    Examples:
        >>> dist = chaospy.Iid(chaospy.Uniform(), 2)
        >>> poly = chaospy.basis(2, 2, dim=2)
        >>> print(poly)
        [q0^2, q0q1, q1^2]
        >>> print(numpy.around(Sens_m_sample(poly, dist, 10000, rule="M"), 4))
        [[0.008  0.0026 0.    ]
         [0.     0.6464 2.1321]]