def Sens_t_sample(poly, dist, samples, rule="R"):
    """
    Total 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 total order sensitivity index for distribution dimensions
            `dim` and polynomial index `pol`.

    Examples:
        >>> dist = chaospy.Iid(chaospy.Uniform(0, 1), 2)
        >>> poly = chaospy.basis(2, 2, dim=2)
        >>> print(poly)
        [q0^2, q0q1, q1^2]
        >>> print(numpy.around(Sens_t_sample(poly, dist, 10000, rule="H"), 4))
        [[ 1.      0.2    -0.3807]
         [ 0.9916  0.9962  1.    ]]
    """
    generator = Saltelli(dist, samples, poly, rule=rule)

    dim = len(dist)
    zeros = [0]*dim
    variance = numpy.var(generator[zeros], -1)
    return numpy.array([
        1-numpy.mean((generator[~index]-generator[zeros])**2, -1,) /
        (2*numpy.where(variance, variance, 1))
        for index in numpy.eye(dim, dtype=bool)
    ])