使用numpy对值进行矢量化计算，而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项可能溢出，从而导致结果无效。你知道吗

我们可以通过使用log空间（使用np.logaddexp求和）来避免这个问题。这个ema_pp_safe比naive方法要贵很多，但仍然比原始循环快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