我的数据是一组N观察到的对,以及它们的频率,即对每对(Xi,Yi),对应于一些Ki,观察到次数(Xi,Yi)。理想情况下,我想计算这些对的所有副本集的Kendall's tau和Spearman's rho,它们由k1+k2+组成+千牛对。问题是k1+k2+…+kn,观察的总数,是巨大的,这样的数据结构将不适合在内存中。
自然地,我考虑了分配第I对的频率,ki/(k1+k2++作为它的权重,计算加权集的秩相关,但我找不到任何工具。在我遇到的秩相关的加权变量(例如,scipy.stats.weightedtau)中,权重表示秩的重要性,而不是对,这与我的原因无关。皮尔逊的r似乎正好有我需要的加权选项,但它不符合我的目的,因为x和y没有线性关系。我想知道我是否遗漏了加权数据点广义相关的一些概念。
到目前为止,我唯一的想法是用一些公共因子c来缩小k1,k2,…,kn,所以第i对拷贝的缩放数量是[ki/c](这里是取整运算符,因为我们需要每对拷贝的整数个数)。通过选择c使得[k1/c]+[k2/c]+…+[kn/c]对可以装入内存,然后我们可以计算出结果集的相关系数tau和rho。然而,ki和kj可以相差很多数量级,因此对于某些ki,c可以显著地大,因此舍入ki/c会导致信息丢失。
upd:可以在指定频率权重的数据集上计算spearman的rho和p值,如下所示:
def frequency_pearsonr(data, frequencies):
"""
Calculates Pearson's r between columns (variables), given the
frequencies of the rows (observations).
:param data: 2-D array with data
:param frequencies: 1-D array with frequencies
:return: 2-D array with pairwise correlations,
2-D array with pairwise p-values
"""
df = frequencies.sum() - 2
Sigma = np.cov(data.T, fweights=frequencies)
sigma_diag = Sigma.diagonal()
Sigma_diag_pairwise_products = np.multiply.outer(sigma_diag, sigma_diag)
# Calculate matrix with pairwise correlations.
R = Sigma / np.sqrt(Sigma_diag_pairwise_products)
# Calculate matrix with pairwise t-statistics. Main diagonal should
# get 1 / 0 = inf.
with np.errstate(divide='ignore'):
T = R / np.sqrt((1 - R * R) / df)
# Calculate matrix with pairwise p-values.
P = 2 * stats.t.sf(np.abs(T), df)
return R, P
def frequency_rank(data, frequencies):
"""
Ranks 1-D data array, given the frequency of each value. Same
values get same "averaged" ranks. Array with ranks is shaped to
match the input data array.
:param data: 1-D array with data
:param frequencies: 1-D array with frequencies
:return: 1-D array with ranks
"""
s = 0
ranks = np.empty_like(data)
# Compute rank for each unique value.
for value in sorted(set(data)):
index_grid = np.ix_(data == value)
# Find total frequency of the value.
frequency = frequencies[index_grid].sum()
ranks[index_grid] = s + 0.5 * (frequency + 1)
s += frequency
return ranks
def frequency_spearmanrho(data, frequencies):
"""
Calculates Spearman's rho between columns (variables), given the
frequencies of the rows (observations).
:param data: 2-D array with data
:param frequencies: 1-D array with frequencies
:return: 2-D array with pairwise correlations,
2-D array with pairwise p-values
"""
# Rank the columns.
ranks = np.empty_like(data)
for i, data_column in enumerate(data.T):
ranks[:, i] = frequency_rank(data_column, frequencies)
# Compute Pearson's r correlation and p-values on the ranks.
return frequency_pearsonr(ranks, frequencies)
# Columns are variables and rows are observations, whose frequencies
# are specified.
data_col1 = np.array([1, 0, 1, 0, 1])
data_col2 = np.array([.67, .25, .75, .2, .6])
data_col3 = np.array([.1, .3, .8, .3, .2])
data = np.array([data_col1, data_col2, data_col3]).T
frequencies = np.array([2, 4, 1, 3, 2])
# Same data, but with observations (rows) actually repeated instead of
# their frequencies being specified.
expanded_data_col1 = np.array([1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1])
expanded_data_col2 = np.array([.67, .67, .25, .25, .25, .25, .75, .2, .2, .2, .6, .6])
expanded_data_col3 = np.array([.1, .1, .3, .3, .3, .3, .8, .3, .3, .3, .2, .2])
expanded_data = np.array([expanded_data_col1, expanded_data_col2, expanded_data_col3]).T
# Compute Spearman's rho for data in both formats, and compare.
frequency_Rho, frequency_P = frequency_spearmanrho(data, frequencies)
Rho, P = stats.spearmanr(expanded_data)
print(frequency_Rho - Rho)
print(frequency_P - P)
上面的特定示例表明,两种方法产生相同的相关性和相同的p值:
[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 1.11022302e-16 0.00000000e+00 -5.55111512e-17]
[ 0.00000000e+00 -5.55111512e-17 0.00000000e+00]]
[[ 0.00000000e+00 -1.35525272e-19 4.16333634e-17]
[ -9.21571847e-19 0.00000000e+00 -5.55111512e-17]
[ 4.16333634e-17 -5.55111512e-17 0.00000000e+00]]
最佳答案
保罗提出的计算肯德尔τ的方法是可行的。不过,您不必将已排序数组的索引指定为列组,未排序数组的索引同样工作良好(如使用加权tau的示例所示)。权重也不需要标准化。
常规(未加权)kendall's tau(在“扩展”数据集上):
stats.kendalltau([0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],
[.25, .25, .25, .25, .2, .2, .2, .667, .667, .75, .6, .6])
KendalltauResult(correlation=0.7977240352174656, pvalue=0.0034446936330652677)
加权kendall's tau(在出现计数为权重的数据集上):
stats.weightedtau([1, 0, 1, 0, 1],
[.667, .25, .75, .2, .6],
rank=False,
weigher=lambda r: [2, 4, 1, 3, 2][r],
additive=False)
WeightedTauResult(correlation=0.7977240352174656, pvalue=nan)
现在,由于weightedtau实现的特殊性,永远不会计算p值。我们可以用最初提供的缩小事件的技巧来近似p值,但我很欣赏其他方法。在我看来,根据可用内存量来决定算法行为似乎很痛苦。
关于python - 在Python中使用频率权重进行排名相关,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/46260215/