python中的numpy数值运算差异与cython加速

Mahone1

各位大佬好，我是一个python小白，最近想把一个数值计算的函数给从pandas版本写成numpy版本，甚至更进阶一点，写成cython化以加速代码运行，代码和我遇到的问题如下：
1、原本的pandas版本的函数是

[Python] 纯文本查看 复制代码

def pandas(data):
    df = pd.DataFrame(data)
    f1 = df.expanding().mean()
    f2 = df.expanding().std()
    shuchu= (df-f1) / f2
    return shuchu

我改写后的numpy版本如下：

[Python] 纯文本查看 复制代码

def numpy(data):
    cumsum = np.cumsum(data)
    index = np.arange(1, len(data) + 1)
    mean_data = np.divide(cumsum, index)
    expanding_std = np.sqrt(np.maximum(0, np.subtract(np.divide(np.cumsum(data**2), index), mean_data**2)))
    shuchu = np.divide(np.subtract(data,mean_data),expanding_std)
    return shuchu

上面这个numpy代码其中计算std标准差的部分在数学公式上是没错的，之所以加了一个和0比较大小的部分是因为一个bug，也就是我想问的一个问题，如下
numpy似乎有精度的问题，0减去0理论上得出的结果应该是0，但是它有时候会计算得出1e-6之类的接近0的极小值，有时候是正数极小值有时候是负数（所以我上面的代码加了一个和0比大小，以免sqrt里的值是负数极小值，会报错）
这个问题，导致数值运算会出错，比如10除一个（0减去0），一般会把计算结果处理成0，但是因为这个精度问题可能就会变成一个大值，比如1e+8之类的
我想问这个问题该怎么解决：
1、我试过用round去压缩精度，但是这会使得其他的数据的计算结果出现微小偏差
2、同时因为我是想为了加速代码的，所以不可能为了精度又去用上一些高精度的数据处理库，这样会大量提高内存和时间的消耗
3、关键是pandas为什么不会出现这种精度问题，比如一个都是0的数据序列去求mean和std，numpy有时候不会出错，有时候会出错，计算出一些极小值，但是pandas比较稳定，目前没看到出错过


总而言之，我想要一个解决方案，解决这个精度问题，同时不要大幅提高内存和时间消耗，请提供一个python代码，或者cython代码，谢谢！


各位大佬快来，我会经常登录吾爱，并且和你讨论细节的，如果你给出的方案是我需要的，我会立刻结算悬赏，跪求各位大佬帮帮忙。

Mahone1

另外，不要完全局限于我给出的代码，只要是能精准计算expanding的mean和std，让原本的数据由此去标准化，就行

chaleaoch

提供一个能重现问题的data例子会帮助到帮助你的人。

helian147

这个AI会帮忙改写

QvQsuipian

在NumPy中处理浮点数时，由于计算机内部表示浮点数的限制，确实可能会遇到精度问题。不过，对于标准差（std）的计算，我们通常可以使用更稳定的算法来避免直接计算(data**2).mean() - mean_data**2导致的精度问题。

在Pandas中，expanding().std()的实现可能使用了更稳定的算法，比如Welford's online algorithm，它可以在不存储整个数据集的情况下计算均值和方差，从而避免了直接计算平方和导致的精度问题。

[Python] 纯文本查看 复制代码

#使用Welford's online algorithm来计算滚动标准差的NumPy版本函数。这个函数应该比你的原始NumPy实现更加稳定和精确
import numpy as np  
  
def numpy_welford(data):  
    n = len(data)  
    mean = np.zeros_like(data, dtype=np.float64)  
    M2 = np.zeros_like(data, dtype=np.float64)  # M2表示方差的无偏估计的累加器  
      
    # 初始化第一个值  
    mean[0] = data[0]  
    M2[0] = 0.0  
      
    # 使用Welford's算法计算滚动均值和方差  
    for i in range(1, n):  
        delta = data[i] - mean[i-1]  
        mean[i] = mean[i-1] + delta / i  
        M2[i] = M2[i-1] + delta * (data[i] - mean[i])  
          
    # 方差的无偏估计需要除以(n-1)，但滚动方差通常除以n  
    std = np.sqrt(M2 / np.arange(1, n + 1, dtype=np.float64))  
      
    # 计算(df - f1) / f2  
    f1 = np.cumsum(data, dtype=np.float64) / np.arange(1, n + 1, dtype=np.float64)  
    shuchu = (data - f1) / std  
      
    # 第一个值的标准差是NaN（或可以设置为0，取决于你的需求）  
    shuchu[0] = np.nan  
      
    return shuchu  
  
# 示例数据  
data = np.random.rand(100)  
result = numpy_welford(data)  
print(result)

Mahone1

QvQsuipian 发表于 2024-5-3 15:24
在NumPy中处理浮点数时，由于计算机内部表示浮点数的限制，确实可能会遇到精度问题。不过，对于标准差（std ...

看起来这个是chatgpt生成的回答，我是不太想用for循环，感觉这样会降低速度，不过我会测试一下你提供的这个代码

QvQsuipian

Mahone1 发表于 2024-5-4 12:57
看起来这个是chatgpt生成的回答，我是不太想用for循环，感觉这样会降低速度，不过我会测试一下你提供的这 ...

数据量有多少？

[Python] 纯文本查看 复制代码

def numpy(data):
    cumsum = np.cumsum(data)
    index = np.arange(1, len(data) + 1)
    mean_data = cumsum / index
    sum_squared = np.cumsum(data**2)
    var_data = np.maximum(0, (sum_squared / index) - mean_data**2)
    std_data = np.sqrt(var_data + np.finfo(float).eps) # 使用eps避免浮点数误差
    shuchu = (data - mean_data) / std_data
    return shuchu

[] 纯文本查看 复制代码

# 保存为numpy_cython.pyx文件

import numpy as np
cimport numpy as np

cpdef numpy_cython(double[:] data):
    cdef int length = data.shape[0]
    cdef double[:] cumsum = np.cumsum(data)
    cdef double[:] mean_data = np.zeros(length)
    cdef double[:] sum_squared = np.cumsum(np.power(data, 2))
    cdef double[:] var_data = np.zeros(length)
    cdef double[:] std_data = np.zeros(length)
    cdef double[:] shuchu = np.zeros(length)
    cdef int i
    
    for i in range(length):
        mean_data[i] = cumsum[i] / (i + 1)
        var_data[i] = max(0, (sum_squared[i] / (i + 1)) - mean_data[i]**2)
        std_data[i] = (var_data[i] + np.finfo(float).eps) ** 0.5
        shuchu[i] = (data[i] - mean_data[i]) / std_data[i]
    
    return shuchu

Mahone1

QvQsuipian 发表于 2024-5-5 00:04
数据量有多少？
[mw_shl_code=python,true]def numpy(data):
    cumsum = np.cumsum(data)

每次传入一维，3000个元素左右的array数据，但是要运算几亿次，所以每一点微小的速度提升对于总体的时间花费而言都是很大的

Mahone1

QvQsuipian 发表于 2024-5-5 00:04
数据量有多少？
[mw_shl_code=python,true]def numpy(data):
    cumsum = np.cumsum(data)

np.finfo(float).eps，这个代码应该对于我的问题没啥用

Mahone1

QvQsuipian 发表于 2024-5-3 15:24
在NumPy中处理浮点数时，由于计算机内部表示浮点数的限制，确实可能会遇到精度问题。不过，对于标准差（std ...

你这个应该是chatgpt的回答，for循环里计算的东西是错误的，比如传入[1,2]，最后一个std居然是0