def _qsturng(p, r, v):
##    print 'q',p
    # r is interpolated through the q to y here we only need to
    # account for when p and/or v are not found in the table.
    global A, p_keys, v_keys

    if p < .1 or p > .999:
        raise ValueError('p must be between .1 and .999')

    if p < .9:
        if v < 2:
            raise ValueError('v must be > 2 when p < .9')
    else:
        if v < 1:
            raise ValueError('v must be > 1 when p >= .9')

    # The easy case. A tabled value is requested.

    #numpy 1.4.1: TypeError: unhashable type: 'numpy.ndarray' :
    p = float(p)
    if isinstance(v, np.ndarray):
        v = v.item()
    if (p,v) in A:
        y = _func(A[(p,v)], p, r, v) + 1.

    elif p not in p_keys and v not in v_keys+([],[1])[p>=.90]:
        # apply bilinear (quadratic) interpolation
        #
        #   p0,v2 +        o         + p1,v2    + p2,v2
        #                    r2
        #
        # 1
        # -                 (p,v)
        # v                x
        #
        #                    r1
        #   p0,v1 +        o         + p1,v1    + p2,v1
        #
        #
        #   p0,v0 +        o r0      + p1,v0    + p2,v0
        #
        #             _ptransform(p)
        #
        # (p1 and v1 may be below or above (p,v). The algorithm
        # works in both cases. For diagramatic simplicity it is
        # shown as above)
        #
        # 1. at v0, v1, and v2 use quadratic interpolation
        #    to find r0, r1, r2
        #
        # 2. use r0, r1, r2 and quadratic interpolaiton
        #    to find y and (p,v)

        # find the 3 closest v values
        v0, v1, v2 = _select_vs(v, p)

        # find the 3 closest p values
        p0, p1, p2 = _select_ps(p)

        # calculate r0, r1, and r2
        r0_sq = _interpolate_p(p, r, v0)**2
        r1_sq = _interpolate_p(p, r, v1)**2
        r2_sq = _interpolate_p(p, r, v2)**2

        # transform v
        v_, v0_, v1_, v2_ = 1./v, 1./v0, 1./v1, 1./v2

        # calculate derivatives for quadratic interpolation
        d2 = 2.*((r2_sq-r1_sq)/(v2_-v1_) - \
                 (r0_sq-r1_sq)/(v0_-v1_)) / (v2_-v0_)
        if (v2_ + v0_) >= (v1_ + v1_):
            d1 = (r2_sq-r1_sq) / (v2_-v1_) - 0.5*d2*(v2_-v1_)
        else:
            d1 = (r1_sq-r0_sq) / (v1_-v0_) + 0.5*d2*(v1_-v0_)
        d0 = r1_sq

        # calculate y
        y = math.sqrt((d2/2.)*(v_-v1_)**2. + d1*(v_-v1_)+ d0)

    elif v not in v_keys+([],[1])[p>=.90]:
        y = _interpolate_v(p, r, v)

    elif p not in p_keys:
        y = _interpolate_p(p, r, v)

    return math.sqrt(2) * -y * \
           scipy.stats.t.isf((1.+p)/2., (v,1e38)[v>1e38])

# make a qsturng functinon that will accept list-like objects
scalar version of qsturng
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/stats/libqsturng/qsturng_.py/_qsturng