How to use the statsmodels.compat.python.range function in statsmodels

def _make_exog_names(exog):
    exog_var = exog.var(0)
    if (exog_var == 0).any():
        # assumes one constant in first or last position
        # avoid exception if more than one constant
        const_idx = exog_var.argmin()
        exog_names = ['x%d' % i for i in range(1, exog.shape[1])]
        exog_names.insert(const_idx, 'const')
    else:
        exog_names = ['x%d' % i for i in range(1, exog.shape[1]+1)]

    return exog_names

Parameters
    ----------
    binsx : array_like, 1d
        binedges

    '''
    binsx = np.asarray(binsx)
    binsy = np.asarray(binsy)
    nx = len(binsx) - 1
    ny = len(binsy) - 1
    probs = np.nan * np.ones((nx, ny)) #np.empty(nx,ny)
    cdf_values = cdf(binsx[:,None], binsy)
    cdf_func = lambda x, y: cdf_values[x,y]
    for xind in range(1, nx+1):
        for yind in range(1, ny+1):
            upper = (xind, yind)
            lower = (xind-1, yind-1)
            #print upper,lower,
            probs[xind-1,yind-1] = prob_bv_rectangle(lower, upper, cdf_func)

    assert not np.isnan(probs).any()
    return probs

##    if sxx is not None and type(sxx) == np.ndarray:
##        sxx_m = sxx[:order+1]
##    else:
##        sxx_m = ut.autocov(s)[:order+1]
    if isacov:
        sxx_m = s
    else:
        sxx_m = acovf(s)[:order+1]  #not tested

    phi = np.zeros((order+1, order+1), 'd')
    sig = np.zeros(order+1)
    # initial points for the recursion
    phi[1,1] = sxx_m[1]/sxx_m[0]
    sig[1] = sxx_m[0] - phi[1,1]*sxx_m[1]
    for k in range(2,order+1):
        phi[k,k] = (sxx_m[k]-np.dot(phi[1:k,k-1], sxx_m[1:k][::-1]))/sig[k-1]
        for j in range(1,k):
            phi[j,k] = phi[j,k-1] - phi[k,k]*phi[k-j,k-1]
        sig[k] = sig[k-1]*(1 - phi[k,k]**2)

    sigma_v = sig[-1]; arcoefs = phi[1:,-1]
    return sigma_v, arcoefs, pacf, phi  #return everything

ibd1 = np.hstack((jj[0:-1][:, None], jj[1:][:, None]))
            ibd1 = [(jj[k], jj[k + 1]) for k in range(len(jj) - 1)]
            ibd.append(ibd1)
        self.ibd = ibd

        # Need to restrict to between-subject pairs
        cpp = []
        for v in model.endog_li:

            # Number of subjects in this group
            m = int(len(v) / self.ncut)

            cpp1 = {}
            # Loop over distinct subject pairs
            for i1 in range(m):
                for i2 in range(i1):
                    # Loop over cut point pairs
                    for k1 in range(self.ncut):
                        for k2 in range(k1+1):
                            if (k2, k1) not in cpp1:
                                cpp1[(k2, k1)] = []
                            j1 = i1*self.ncut + k1
                            j2 = i2*self.ncut + k2
                            cpp1[(k2, k1)].append([j2, j1])

            for k in cpp1.keys():
                cpp1[k] = np.asarray(cpp1[k])
            cpp.append(cpp1)

        self.cpp = cpp

        # Initialize the dependence parameters

def edgeworth(intervals):
            """Compute the Edgeworth expansion term of Sison & Glaz's formula
            (1) (approximated probability for multinomial proportions in a
            given box)."""
            # Compute means and central moments of the truncated poisson
            # variables.
            mu_r1, mu_r2, mu_r3, mu_r4 = [
                np.array([truncated_poisson_factorial_moment(interval, r, p)
                          for (interval, p) in zip(intervals, counts)])
                for r in range(1, 5)
            ]
            mu = mu_r1
            mu2 = mu_r2 + mu - mu ** 2
            mu3 = mu_r3 + mu_r2 * (3 - 3 * mu) + mu - 3 * mu ** 2 + 2 * mu ** 3
            mu4 = (mu_r4 + mu_r3 * (6 - 4 * mu) +
                   mu_r2 * (7 - 12 * mu + 6 * mu ** 2) +
                   mu - 4 * mu ** 2 + 6 * mu ** 3 - 3 * mu ** 4)

            # Compute expansion factors, gamma_1 and gamma_2.
            g1 = mu3.sum() / mu2.sum() ** 1.5
            g2 = (mu4.sum() - 3 * (mu2 ** 2).sum()) / mu2.sum() ** 2

            # Compute the expansion itself.
            x = (n - mu.sum()) / np.sqrt(mu2.sum())
            phi = np.exp(- x ** 2 / 2) / np.sqrt(2 * np.pi)
            H3 = x ** 3 - 3 * x

A = copy.copy(A)
        A.namespace = namespace
        B = copy.copy(B)
        B.namespace = namespace

    a = A(values=True)[0]
    b = B(values=True)[0]

    if len(a) != len(b):
        raise ValueError('A() and B() should be sequences of the same length')

    nA = len(set(a))
    nB = len(set(b))
    n = max(nA, nB)

    AB = [(a[i],b[i]) for i in range(len(a))]
    nAB = len(set(AB))

    if nAB == n:
        if nA > nB:
            F = A
        else:
            F = B
        return (True, F)
    else:
        return (False, None)

alpha = np.exp(params[-1])
        else:
            alpha = params[-1]
        a1 = 1/alpha
        params = params[:-1]

        exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))
        const_arr = a1*mu*(a1+y)/(mu+a1)**2
        for i in range(dim):
            for j in range(dim):
                if j > i:
                    continue
                hess_arr[i,j] = np.sum(-exog[:,i,None] * exog[:,j,None] *
                                       const_arr, axis=0)
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        # for dl/dparams dalpha
        da1 = -alpha**-2
        dldpda = np.sum(mu*exog*(y-mu)*da1/(mu+a1)**2 , axis=0)
        hess_arr[-1,:-1] = dldpda
        hess_arr[:-1,-1] = dldpda

        # for dl/dalpha dalpha
        #NOTE: polygamma(1,x) is the trigamma function
        da2 = 2*alpha**-3

exog = self.exog
        y = self.endog[:,None]
        mu = self.predict(params)[:,None]
        mu_p = np.power(mu, p)
        a1 = 1 + alpha * mu_p
        a2 = mu + alpha * mu_p * y
        a3 = alpha * p * mu ** (p - 1)
        a4 = a3 * y
        a5 = p * mu ** (p - 1)
        dmudb = mu * exog

        # for dl/dparams dparams
        dim = exog.shape[1]
        hess_arr = np.empty((dim+1,dim+1))

        for i in range(dim):
            for j in range(i + 1):
                hess_arr[i,j] = np.sum(mu * exog[:,i,None] * exog[:,j,None] *
                    (mu * (a3 * a4 / a1**2 - 2 * a3**2 * a2 / a1**3 + 2 * a3 *
                    (a4 + 1) / a1**2 - a4 * p / (mu * a1) + a3 * p * a2 /
                    (mu * a1**2) + a4 / (mu * a1) - a3 * a2 / (mu * a1**2) +
                    (y - 1) * a4 * (p - 1) / (a2 * mu) - (y - 1) *
                    (1 + a4)**2 / a2**2 - a4 * (p - 1) / (a1 * mu) - 1 /
                    mu**2) + (-a4 / a1 + a3 * a2 / a1**2 + (y - 1) *
                    (1 + a4) / a2 - (1 + a4) / a1 + 1 / mu)), axis=0)
        tri_idx = np.triu_indices(dim, k=1)
        hess_arr[tri_idx] = hess_arr.T[tri_idx]

        # for dl/dparams dalpha
        dldpda = np.sum((2 * a4 * mu_p / a1**2 - 2 * a3 * mu_p * a2 / a1**3 -
                        mu_p * y * (y - 1) * (1 + a4) / a2**2 + mu_p *
                        (1 + a4) / a1**2 + a5 * y * (y - 1) / a2 - 2 *

"""

        endog = self.model.endog_li
        cpp = self.cpp
        cached_means = self.model.cached_means

        # This will happen if all the clusters have only
        # one observation
        if len(cpp[0]) == 0:
            return

        tables = {}
        for ii in cpp[0]:
            tables[ii] = np.zeros((2, 2), dtype=np.float64)

        for i in range(self.model.num_group):

            endog_expval, _ = cached_means[i]

            emat_11 = self.get_eyy(endog_expval, i)
            emat_10 = endog_expval[:, None] - emat_11
            emat_01 = -emat_11 + endog_expval
            emat_00 = 1. - (emat_11 + emat_10 + emat_01)

            cpp1 = cpp[i]
            for ky in iterkeys(cpp1):
                ix = cpp1[ky]
                tables[ky][1, 1] += emat_11[ix[:, 0], ix[:, 1]].sum()
                tables[ky][1, 0] += emat_10[ix[:, 0], ix[:, 1]].sum()
                tables[ky][0, 1] += emat_01[ix[:, 0], ix[:, 1]].sum()
                tables[ky][0, 0] += emat_00[ix[:, 0], ix[:, 1]].sum()

old version without precomputing cdf values

    Parameters
    ----------
    binsx : array_like, 1d
        binedges

    '''
    binsx = np.asarray(binsx)
    binsy = np.asarray(binsy)
    nx = len(binsx) - 1
    ny = len(binsy) - 1
    probs = np.nan * np.ones((nx, ny)) #np.empty(nx,ny)
    for xind in range(1, nx+1):
        for yind in range(1, ny+1):
            upper = (binsx[xind], binsy[yind])
            lower = (binsx[xind-1], binsy[yind-1])
            #print upper,lower,
            probs[xind-1,yind-1] = prob_bv_rectangle(lower, upper, cdf)

    assert not np.isnan(probs).any()
    return probs