Скользящее сравнение между значением и прошедшим окном с процентилем / квантилем

Question

Я бы хотел сравнить каждое значение x массива со скользящим окном из n предыдущих значений. Точнее, я хотел бы видеть , в каком процентиле было бы это новое значение x, если бы мы добавили его в предыдущее окно :

import numpy as np
A = np.array([1, 4, 9, 28, 28.5, 2, 283, 3.2, 7, 15])
print A
n = 4  # window width
for i in range(len(A)-n):
    W = A[i:i+n]
    x = A[i+n]
    q = sum(W <= x) * 1.0 / n
    print 'Value:', x, ' Window before this value:', W, ' Quantile:', q

[1. 4. 9. 28. 28.5 2. 283. 3.2 7. 15.]
Значение: 28,5. Окно до этого значения: [1. 4. 9. 28.] Квантиль: 1.0
Значение: 2.0 Окно до этого значения: [4. 9. 28. 28.5] Квантиль: 0.0
Значение: 283,0. Окно до этого значения: [9. 28. 28,5 2.] Квантиль: 1,0
Значение: 3.2 Окно до этого значения: [28. 28.5 2. 283.] Квантиль: 0.25
Значение: 7.0. Окно до этого значения: [28,5 2. 283. 3.2] Квантиль: 0,5
Значение: 15.0. Окно до этого значения: [2. 283. 3.2 7.] Квантиль: 0.75

Вопрос: Как называется это вычисление? Есть ли хитрый способ вычислить это более эффективно для массивов миллионов элементов (с n, которое может быть ~ 5000)?

---

Примечание: вот симуляция для 1М предметов и n = 5000, но это займет ~ 2 часа:

import numpy as np
A = np.random.random(1000*1000)  # the following is not very interesting with a [0,1]
n = 5000                         # uniform random variable, but anyway...
Q = np.zeros(len(A)-n)
for i in range(len(Q)):
    Q[i] = sum(A[i:i+n] <= A[i+n]) * 1.0 / n
    if i % 100 == 0: 
        print "%.2f %% already done. " % (i * 100.0 / len(A))

print Q

Примечание: это не похоже на Как вычислить движущийся (или вращающийся, если хотите) процентиль / квантиль для массива 1d в numpy?

Answer 1

Дополнительный тест: сравнение между этим решением и этим решением :

import numpy as np, time

A = np.random.random(1000*1000)
n = 5000

def compare_strides (arr, n):
   return (np.lib.stride_tricks.as_strided(arr, shape=(n,arr.size-n), strides=(arr.itemsize,arr.itemsize)) <= arr[n:]).sum(0)

# Test #1: with strides ===> 11.0 seconds
t0 = time.time()
nb_chunk = 10*1000
Q = np.concatenate([compare_strides(A[chunk*nb_chunk:(chunk+1)*nb_chunk+n],n) for chunk in range(0,A[n:].size/nb_chunk+1)])/(1.*n)
print time.time() - t0, Q

# Test #2: with just np.sum ===> 18.0 seconds
t0 = time.time()
Q2 = np.zeros(len(A)-n)
for i in range(len(Q2)):
    Q2[i] = np.sum(A[i:i+n] <= A[i+n])
Q2 *= 1.0 / n  # here the multiplication is vectorized; if instead, we move this multiplication to the previous line: np.sum(A[i:i+n] <= A[i+n]) * 1.0 / n, it is 6 seconds slower
print time.time() - t0, Q2

print all(Q == Q2)

Есть и другой (лучший) способ с декораторами numba и @jit. Тогда это намного быстрее: всего 5,4 секунды !

from numba import jit
import numpy as np

@jit  # if you remove this line, it is much slower (similar to Test #2 above)
def doit():
    A = np.random.random(1000*1000)
    n = 5000
    Q2 = np.zeros(len(A)-n)
    for i in range(len(Q2)):
        Q2[i] = np.sum(A[i:i+n] <= A[i+n])
    Q2 *= 1.0/n
    print(Q2)

doit()

При добавлении распараллеливания numba это происходит еще быстрее: 1.8 секунд!

import numpy as np
from numba import jit, prange

@jit(parallel=True)
def doit(A, Q, n):
    for i in prange(len(Q)):
        Q[i] = np.sum(A[i:i+n] <= A[i+n])

A = np.random.random(1000*1000)
n = 5000
Q = np.zeros(len(A)-n)    
doit(A, Q, n)

Answer 2

Использование np.sum вместо суммы уже упоминалось, поэтому мое единственное оставленное предложение - дополнительно рассмотреть возможность использования pandas и функции скользящего окна, к которой вы можете применить любую произвольную функцию:

import numpy as np
import pandas as pd

A = np.random.random(1000*1000)
df = pd.DataFrame(A)
n = 5000

def fct(x):
    return np.sum(x[:-1] <= x[-1]) * 1.0 / (len(x)-1)

percentiles = df.rolling(n+1).apply(fct)
print(percentiles)

Answer 3

вы можете использовать np.lib.stride_tricks.as_strided, как в принятом ответе на вопрос, который вы связали. С первым примером, который вы даете, это довольно легко понять:

A = np.array([1, 4, 9, 28, 28.5, 2, 283, 3.2, 7, 15])
n=4
print (np.lib.stride_tricks.as_strided(A, shape=(A.size-n,n),
                                       strides=(A.itemsize,A.itemsize)))
# you get the A.size-n columns of the n rolling elements
array([[  1. ,   4. ,   9. ,  28. ,  28.5,   2. ],
       [  4. ,   9. ,  28. ,  28.5,   2. , 283. ],
       [  9. ,  28. ,  28.5,   2. , 283. ,   3.2],
       [ 28. ,  28.5,   2. , 283. ,   3.2,   7. ]])

Теперь, чтобы выполнить вычисление, вы можете сравнить этот массив с A [n:], sum по строкам и разделить на n:

print ((np.lib.stride_tricks.as_strided(A, shape=(n,A.size-n),
                                        strides=(A.itemsize,A.itemsize)) 
          <= A[n:]).sum(0)/(1.*n))
[1.   0.   1.   0.25 0.5  0.75] # same anwser

Теперь проблема заключается в размере ваших данных (несколько M и n около 5000), но вы не уверены, что можете использовать этот метод напрямую. Одним из способов может быть разделение данных на части. Давайте определим функцию

def compare_strides (arr, n):
   return (np.lib.stride_tricks.as_strided(arr, shape=(n,arr.size-n),
                                           strides=(arr.itemsize,arr.itemsize)) 
            <= arr[n:]).sum(0)

и сделайте кусок, с np.concatenate и не забудьте разделить на n:

nb_chunk = 1000 #this number depends on the capacity of you computer, 
                # not sure how to optimize it
Q = np.concatenate([compare_strides(A[chunk*nb_chunk:(chunk+1)*nb_chunk+n],n) 
                    for chunk in range(0,A[n:].size/nb_chunk+1)])/(1.*n)

Я не могу выполнить тест 1M - 5000, но на 5000 - 100 вижу разницу в timeit:

A = np.random.random(5000)
n = 100

%%timeit
Q = np.zeros(len(A)-n)
for i in range(len(Q)):
    Q[i] = sum(A[i:i+n] <= A[i+n]) * 1.0 / n

#1 loop, best of 3: 6.75 s per loop

%%timeit
nb_chunk = 100
Q1 = np.concatenate([compare_strides(A[chunk*nb_chunk:(chunk+1)*nb_chunk+n],n) 
                    for chunk in range(0,A[n:].size/nb_chunk+1)])/(1.*n)

#100 loops, best of 3: 7.84 ms per loop

#check for egality
print ((Q == Q1).all())
Out[33]: True

См. Разницу во времени от 6750 мс до 7,84 мс. Надеюсь, что это работает на больших данных

Answer 4

Ваш код очень медленный, потому что вы используете собственный Python sum() вместо numpy.sum() или numpy.array.sum(); sum() Python должен преобразовать все необработанные значения в объекты Python перед выполнением вычислений, что действительно медленно. Просто изменив sum(...) на np.sum(...) или (...).sum(), время выполнения упадет до 20 секунд.

Answer 5

Вы можете использовать np.quantile вместо sum(A[i:i+n] <= A[i+n]) * 1.0 / n. Это может быть так же хорошо, как и получается. Не уверен, что действительно есть лучший подход к вашему вопросу.