regionsRanker.py 文件源码

def lleval(self, Uy, UX, Sd, yKy, logdetK, logdetXX, logdelta, UUXUUX, UUXUUy, UUyUUy, numIndividuals, reml):
        N = numIndividuals
        D = UX.shape[1]

        UXS = UX / np.lib.stride_tricks.as_strided(Sd, (Sd.size, D), (Sd.itemsize,0))
        XKy = UXS.T.dot(Uy)         
        XKX = UXS.T.dot(UX) 

        if (Sd.shape[0] < numIndividuals):
            delta = np.exp(logdelta)
            denom = delta
            XKX += UUXUUX / denom
            XKy += UUXUUy / denom
            yKy += UUyUUy / denom           
            logdetK += (numIndividuals-Sd.shape[0]) * logdelta      

        [SxKx,UxKx]= la.eigh(XKX)   
        i_pos = SxKx>1E-10
        beta = np.dot(UxKx[:,i_pos], (np.dot(UxKx[:,i_pos].T, XKy) / SxKx[i_pos]))
        r2 = yKy-XKy.dot(beta)

        if reml:
            logdetXKX = np.log(SxKx).sum()
            sigma2 = (r2 / (N - D))
            ll =  -0.5 * (logdetK + (N-D)*np.log(2.0*np.pi*sigma2) + (N-D) + logdetXKX - logdetXX)
        else:
            sigma2 = r2 / N
            ll =  -0.5 * (logdetK + N*np.log(2.0*np.pi*sigma2) + N)

        return ll, sigma2, beta, r2