使用 numpy 对值进行矢量化计算，这需要先前计算的值

Question

我正在尝试从 Investopedia 计算一个特定的 EMA 公式，看起来像

EmaToday = (ValueToday ∗ (Smoothing / 1+Days)) 
           + (EmaYesterday * (1 - (Smoothing / 1+Days)))

我们可以将其简化为：

Smoothing and Days are constants.
Let's call (Smoothing / 1 + Days) as 'M'

The simplified equation becomes:
EmaToday = ((ValueToday - EmaYesterday) * M) + EmaYesterday

我们可以在传统的 python 中使用如下循环来做到这一点：

# Initialize an empty numpy array to hold calculated ema values
emaTodayArray = np.zeros((1, valueTodayArray.size - Days), dtype=np.float32)

ema = emaYesterday
# Calculate ema
for i, valueToday in enumerate(np.nditer(valueList)):
    ema = ((valueToday - ema) * M) + ema
    emaTodayArray[i] = ema

emaTodayArray 保存所有计算的 EMA 值。

我很难弄清楚如何将其完全矢量化，因为每次新计算都需要 emaYesterday 值。

如果首先可以使用 numpy 进行完全矢量化，如果有人能告诉我方法，我将不胜感激。

​

Answer 1

注意：我必须填写一些假人才能使您的代码运行，请检查它们是否正常。

可以通过转换 ema[i] ~> ema'[i] = ema[i] x (1-M)^-i 来对循环进行矢量化，然后它就变成了 cumsum。

这在下面被实现为ema_pp_naive。

这种方法的问题在于，对于中等大小的i (~10^3)​​，(1-M)^-i 项可能会溢出，导致结果无用。

我们可以通过转到日志空间来规避这个问题（使用np.logaddexp 进行求和）。这个ema_pp_safe 比天真的方法要贵很多，但仍然比原始循环快 10 倍以上。在我快速而肮脏的测试中，这给出了一百万个甚至更多的正确结果。

代码：

import numpy as np

K = 1000
Days = 0

emaYesterday = np.random.random()
valueTodayArray = np.random.random(K)
M = np.random.random()

valueList = valueTodayArray


import time

T = []

T.append(time.perf_counter())

# Initialize an empty numpy array to hold calculated ema values
emaTodayArray = np.zeros((valueTodayArray.size - Days), dtype=np.float32)

ema = emaYesterday
# Calculate ema
for i, valueToday in enumerate(np.nditer(valueList)):
    ema = ((valueToday - ema) * M) + ema
    emaTodayArray[i] = ema

T.append(time.perf_counter())

scaling = np.broadcast_to(1/(1-M), valueTodayArray.size+1).cumprod()
ema_pp_naive = ((np.concatenate([[emaYesterday], valueTodayArray * M]) * scaling).cumsum() / scaling)[1:]

T.append(time.perf_counter())

logscaling = np.log(1-M)*np.arange(valueTodayArray.size+1)
log_ema_pp = np.logaddexp.accumulate(np.log(np.concatenate([[emaYesterday], valueTodayArray * M])) - logscaling) + logscaling
ema_pp_safe = np.exp(log_ema_pp[1:])

T.append(time.perf_counter())

print(f'K = {K}')
print('naive method correct:', np.allclose(ema_pp_naive, emaTodayArray))
print('safe method correct:', np.allclose(ema_pp_safe, emaTodayArray))
print('OP {:.3f} ms   naive {:.3f} ms   safe {:.3f} ms'.format(*np.diff(T)*1000))

示例运行：

K = 100
naive method correct: True
safe method correct: True
OP 0.236 ms   naive 0.061 ms   safe 0.053 ms

K = 1000
naive method correct: False
safe method correct: True
OP 2.397 ms   naive 0.224 ms   safe 0.183 ms

K = 1000000
naive method correct: False
safe method correct: True
OP 2145.956 ms   naive 18.342 ms   safe 108.528 ms