将两组相同的数据传递给自制 pearson 的相关实现时返回 0.999...2
Homemade pearson's correlation implementation returning 0.999...2 when passing two of the same sets of data to it
我受够了 scipy 和 numpy,并决定继续并根据某个地方的 SO 答案继续进行另一个实现。
from statistics import pstdev, mean
def pearson(x, y):
sx = []
sy = []
mx = mean(x)
my = mean(y)
stdx = pstdev(x)
stdy = pstdev(y)
for i in x:
sx.append((i - mx) / stdx)
for j in y:
sy.append((j - my) / stdy)
return sum([i * j for i, j in zip(sx, sy)]) / len(x)
我给它传了几个数字,看它给的是不是和scipy.stats.pearsonr
一样的东西,好像没问题。最后的一些数字有所不同,但没有什么开创性的......
直到我尝试将同一组数据传递给它作为 x
和 y
。当我这样做时,我返回 0.9999999999999992
,当 scipy 和 numpy 都说它是 1.0
.
这个实现有问题吗?我正在使用 population stdev 而不是样本 stdev,据我所知,numpy 和 scipy 都使用它。我想知道为什么它没有按应有的方式返回 1.0
。可能是 python 本身的浮动问题吗?我已经在 Py 2 和 3 中尝试过,并且我在两者中都得到了 0.999...
。
如果需要,我传入的数据集是:
[7, 1, 5, 1, 8, 5, 9, 8, 5, 10, 5, 8, 1, 8, 8, 8, 10, 4, 8, 9, 9, 6, 8, 7, 8, 5, 10, 5, 6, 8, 8, 7, 9, 4, 6, 10, 7, 10, 4, 5, 4, 7, 4, 8, 9, 10, 9, 8, 7, 8, 6, 8, 6, 6, 5, 7, 7, 7, 7, 3, 7, 8, 6, 8, 5, 7, 8, 7, 8, 6, 8, 6, 9, 6, 6, 6, 8, 9, 5, 7, 9, 2, 9, 6, 7, 6, 7, 7, 5, 5, 7, 7, 8, 6, 9, 1, 3, 6, 7, 9, 7, 7, 6, 9, 9, 4, 9, 9, 7, 9, 6, 2, 2, 8, 4, 7, 7, 6, 3, 7, 3, 5, 10, 9, 8, 10, 8, 7, 4, 7, 8, 9, 8, 4, 7, 9, 7, 7, 6, 8, 8, 9, 9, 7, 4, 4, 7, 3, 9, 3, 1, 8, 3, 9, 4, 8, 3, 9, 8, 8, 7, 9, 9, 8, 10, 8, 3, 10, 4, 7, 7, 10, 8, 7, 8, 7, 1, 8, 9, 5, 7, 5, 5, 3, 5, 7, 7, 7, 2, 4, 1, 6, 9, 9, 7, 7, 10, 9, 2, 9, 8, 2, 5, 1, 2, 5, 9, 1, 4, 8, 9, 6, 4, 4, 7, 3, 7, 9, 4, 3, 7, 8, 7, 6, 8, 8, 7]
您对浮点行为的期望过于乐观。根据经验,您不会对结果不完全是 1.0 感到惊讶。例如,试试这个小得多的输入:
[7, 1, 5]
在我的盒子上,你的函数 returns 1.0000000000000002。 "Close to" 1.0,但不完全是 1.0。一般来说,这是您所能期望的最好结果。
想知道为什么,想想这个 "should" 计算的是什么:
math.sqrt(x)**2 == x
在数学上(以无限精度计算),这应该总是 return 正确。但是在浮点数中(无论使用多少精度,只要精度有界),它不可能总是为真。事实上,反例很容易找到;比如,刚才在我的盒子上:
>>> math.sqrt(2)**2
2.0000000000000004
问题在于,具有有限精度,sqrt()
必然是多对一函数。它将域 1..N
压缩到范围 1..sqrt(N)
中,并且在有限精度下,域的基数大于范围的基数。因此,域中必须存在不同的 x
和 y
映射到范围内的相同值,因此不存在精确的函数逆。
你的函数比普通函数更复杂sqrt
,但相同的原理在起作用。
是的,这与浮点行为有关。您可以尝试使用 decimal 模块
from decimal import *
data = [7, 1, 5, 1, 8, 5, 9, 8, 5, 10, 5, 8, 1, 8, 8, 8, 10, 4, 8]
data = [Decimal(x) for x in data]
print(pearson(data, data))
请注意,您还需要使用小数来计算平均值和标准差。 您可以使用 x ** Decimal('0.5')
而不是 sqrt(x)
使用 decimal 的 sqrt 函数,正如 Tim Peters 在评论中所解释的那样。
我受够了 scipy 和 numpy,并决定继续并根据某个地方的 SO 答案继续进行另一个实现。
from statistics import pstdev, mean
def pearson(x, y):
sx = []
sy = []
mx = mean(x)
my = mean(y)
stdx = pstdev(x)
stdy = pstdev(y)
for i in x:
sx.append((i - mx) / stdx)
for j in y:
sy.append((j - my) / stdy)
return sum([i * j for i, j in zip(sx, sy)]) / len(x)
我给它传了几个数字,看它给的是不是和scipy.stats.pearsonr
一样的东西,好像没问题。最后的一些数字有所不同,但没有什么开创性的......
直到我尝试将同一组数据传递给它作为 x
和 y
。当我这样做时,我返回 0.9999999999999992
,当 scipy 和 numpy 都说它是 1.0
.
这个实现有问题吗?我正在使用 population stdev 而不是样本 stdev,据我所知,numpy 和 scipy 都使用它。我想知道为什么它没有按应有的方式返回 1.0
。可能是 python 本身的浮动问题吗?我已经在 Py 2 和 3 中尝试过,并且我在两者中都得到了 0.999...
。
如果需要,我传入的数据集是:
[7, 1, 5, 1, 8, 5, 9, 8, 5, 10, 5, 8, 1, 8, 8, 8, 10, 4, 8, 9, 9, 6, 8, 7, 8, 5, 10, 5, 6, 8, 8, 7, 9, 4, 6, 10, 7, 10, 4, 5, 4, 7, 4, 8, 9, 10, 9, 8, 7, 8, 6, 8, 6, 6, 5, 7, 7, 7, 7, 3, 7, 8, 6, 8, 5, 7, 8, 7, 8, 6, 8, 6, 9, 6, 6, 6, 8, 9, 5, 7, 9, 2, 9, 6, 7, 6, 7, 7, 5, 5, 7, 7, 8, 6, 9, 1, 3, 6, 7, 9, 7, 7, 6, 9, 9, 4, 9, 9, 7, 9, 6, 2, 2, 8, 4, 7, 7, 6, 3, 7, 3, 5, 10, 9, 8, 10, 8, 7, 4, 7, 8, 9, 8, 4, 7, 9, 7, 7, 6, 8, 8, 9, 9, 7, 4, 4, 7, 3, 9, 3, 1, 8, 3, 9, 4, 8, 3, 9, 8, 8, 7, 9, 9, 8, 10, 8, 3, 10, 4, 7, 7, 10, 8, 7, 8, 7, 1, 8, 9, 5, 7, 5, 5, 3, 5, 7, 7, 7, 2, 4, 1, 6, 9, 9, 7, 7, 10, 9, 2, 9, 8, 2, 5, 1, 2, 5, 9, 1, 4, 8, 9, 6, 4, 4, 7, 3, 7, 9, 4, 3, 7, 8, 7, 6, 8, 8, 7]
您对浮点行为的期望过于乐观。根据经验,您不会对结果不完全是 1.0 感到惊讶。例如,试试这个小得多的输入:
[7, 1, 5]
在我的盒子上,你的函数 returns 1.0000000000000002。 "Close to" 1.0,但不完全是 1.0。一般来说,这是您所能期望的最好结果。
想知道为什么,想想这个 "should" 计算的是什么:
math.sqrt(x)**2 == x
在数学上(以无限精度计算),这应该总是 return 正确。但是在浮点数中(无论使用多少精度,只要精度有界),它不可能总是为真。事实上,反例很容易找到;比如,刚才在我的盒子上:
>>> math.sqrt(2)**2
2.0000000000000004
问题在于,具有有限精度,sqrt()
必然是多对一函数。它将域 1..N
压缩到范围 1..sqrt(N)
中,并且在有限精度下,域的基数大于范围的基数。因此,域中必须存在不同的 x
和 y
映射到范围内的相同值,因此不存在精确的函数逆。
你的函数比普通函数更复杂sqrt
,但相同的原理在起作用。
是的,这与浮点行为有关。您可以尝试使用 decimal 模块
from decimal import *
data = [7, 1, 5, 1, 8, 5, 9, 8, 5, 10, 5, 8, 1, 8, 8, 8, 10, 4, 8]
data = [Decimal(x) for x in data]
print(pearson(data, data))
请注意,您还需要使用小数来计算平均值和标准差。 您可以使用 使用 decimal 的 sqrt 函数,正如 Tim Peters 在评论中所解释的那样。x ** Decimal('0.5')
而不是 sqrt(x)