ime by providing a NumPy array or list for `df`:

    >>> stdtrit([1, 2, 3], 0.7)
    array([0.72654253, 0.6172134 , 0.58438973])

    It is possible to calculate the function at several points for several
    different degrees of freedom simultaneously by providing arrays for `df`
    and `p` with shapes compatible for broadcasting. Compute `stdtrit` at
    4 points for 3 degrees of freedom resulting in an array of shape 3x4.

    >>> dfs = np.array([[1], [2], [3]])
    >>> p = np.array([0.2, 0.4, 0.7, 0.8])
    >>> dfs.shape, p.shape
    ((3, 1), (4,))

    >>> stdtrit(dfs, p)
    array([[-1.37638192, -0.3249197 ,  0.72654253,  1.37638192],
           [-1.06066017, -0.28867513,  0.6172134 ,  1.06066017],
           [-0.97847231, -0.27667066,  0.58438973,  0.97847231]])

    The t distribution is also available as `scipy.stats.t`. Calling `stdtrit`
    directly can be much faster than calling the ``ppf`` method of
    `scipy.stats.t`. To get the same results, one must use the following
    parametrization: ``scipy.stats.t(df).ppf(x) = stdtrit(df, x)``.

    >>> from scipy.stats import t
    >>> df, x = 3, 0.5
    >>> stdtrit_result = stdtrit(df, x)  # this can be faster than below
    >>> stats_result = t(df).ppf(x)
    >>> stats_result == stdtrit_result  # test that results are equal
    True
    Ú