smtpred.py 文件源码

def ols_weights(n, h2, rg, mtot=1e6):
    """
    Args: n: vector with sample size for each trait
          h2: vector with SNP heritabilities
          rg: vector with rg for each pair of traits (3 traits: 1,2; 1,3; 2,3)
          mtot: total number of markers (doesn't change result much)
    Returns: ntraits * ntraits array with ols weights. weights in each row are for are for a multi-trait predictor of the trait in this row
    """
    ntraits = len(n)
    gcovmat = get_gcovmat(h2, rg)
    print(gcovmat)
    V = gcovmat / mtot
    numpy.fill_diagonal(V, ols_variances(n, h2, mtot))
    C = gcovmat / mtot

    weights = numpy.zeros([ntraits, ntraits])
    for i in range(ntraits):
        nonzero = V[i,] != 0
        Vi = V[numpy.array(numpy.where(nonzero)[0])[:, None], nonzero]
        Vinv = numpy.linalg.inv(Vi)
        weights[i, nonzero] = numpy.dot(Vinv, C[i, nonzero])
    print(weights)
    return weights