[python] NumPy를 사용하여 이동 평균을 계산하는 방법은 무엇입니까?

단순히 numpy / scipy의 이동 평균을 계산하는 함수가 없어 복잡한 솔루션으로 이어지는 것 같습니다 .

내 질문은 두 가지입니다.

  • numpy로 이동 평균을 (올바르게) 구현하는 가장 쉬운 방법은 무엇입니까?
  • 이것은 사소하지 않고 오류가 발생하기 쉬운 것처럼 보이므로이 경우에 배터리를 포함 하지 않는 것이 좋은 이유가 있습니까?


답변

그냥 간단합니다 이동 평균 비를 가중, 당신은 쉽게 그것을 구현할 수 np.cumsum있는, 할 수는 있습니다 빠른 FFT에 비해 기반 방법 :

EDIT 코드에서 Bean이 발견 한 잘못된 인덱싱을 하나씩 수정했습니다. 편집하다

def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

>>> a = np.arange(20)
>>> moving_average(a)
array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
        12.,  13.,  14.,  15.,  16.,  17.,  18.])
>>> moving_average(a, n=4)
array([  1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
        10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5])

그래서 대답은 구현하기가 정말 쉽고 numpy는 이미 특수한 기능으로 약간 부풀어 오른 것입니다.


답변

NumPy의 특정 도메인 특정 기능이 부족한 것은 아마도 코어 팀의 원칙과 NumPy의 주요 지시문에 대한 충실 성 때문일 것입니다 . N- 차원 배열 유형을 제공하고 해당 배열을 만들고 색인화하는 기능을 제공합니다. 많은 기본 목표와 마찬가지로 이것은 작지 않으며 NumPy는 훌륭하게 수행합니다.

(훨씬) 더 큰 SciPy 에는 훨씬 더 큰 도메인 별 라이브러리 모음 ( SciPy 개발자의 하위 패키지 라고 함 )이 포함되어 있습니다 ( 예 : 수치 최적화 ( optimize ), 신호 처리 ( signal ) 및 적분 계산 ( 통합 )).

내 생각에 당신이 추구 하는 기능은 적어도 하나의 SciPy 서브 패키지 ( scipy.signal 아마도); 그러나 나는 먼저 SciPy scikits 컬렉션을 살펴보고 관련 scikit (s)를 식별하고 거기에서 관심있는 기능을 찾습니다.

Scikits 는 NumPy / SciPy를 기반으로 독립적으로 개발 된 패키지이며 특정 기술 분야 (예 : scikits-image , scikits-learn 등)를 대상으로합니다. 이들 중 일부는 (특히 수치 최적화를위한 멋진 OpenOpt ) 높은 평가를 받았습니다. 비교적 새로운 scikits 루 브릭 하에 거주하기로 선택하기 훨씬 전에 성숙한 프로젝트 . Scikits은 (30) 등에 대한 위의 목록에 좋아 한 적이 홈페이지 scikits 적어도 그 중 몇 가지가 더 이상 활성 개발중인 비록.

이 조언을 따르면 scikits-timeseries로 연결됩니다 . 그러나 해당 패키지는 더 이상 활발하게 개발되지 않습니다. 실제로 Pandas사실상 NumPy 기반 시계열 라이브러리 인 AFAIK가되었습니다 .

Pandas 에는 이동 평균 을 계산하는 데 사용할 수있는 여러 함수가 있습니다 . 이들 중 가장 간단한 것은 아마도 rolling_mean 일 것입니다 .

>>> # the recommended syntax to import pandas
>>> import pandas as PD
>>> import numpy as NP

>>> # prepare some fake data:
>>> # the date-time indices:
>>> t = PD.date_range('1/1/2010', '12/31/2012', freq='D')

>>> # the data:
>>> x = NP.arange(0, t.shape[0])

>>> # combine the data & index into a Pandas 'Series' object
>>> D = PD.Series(x, t)

이제 Series 객체와 창 크기를 전달하는 rolling_mean 함수를 호출하십시오 . 아래 예제에서는 10 일 입니다.

>>> d_mva = PD.rolling_mean(D, 10)

>>> # d_mva is the same size as the original Series
>>> d_mva.shape
    (1096,)

>>> # though obviously the first w values are NaN where w is the window size
>>> d_mva[:3]
    2010-01-01         NaN
    2010-01-02         NaN
    2010-01-03         NaN

효과가 있는지 확인합니다. 예를 들어 원래 시리즈의 값 10-15를 롤링 평균으로 평활화 된 새 시리즈와 비교

>>> D[10:15]
     2010-01-11    2.041076
     2010-01-12    2.041076
     2010-01-13    2.720585
     2010-01-14    2.720585
     2010-01-15    3.656987
     Freq: D

>>> d_mva[10:20]
      2010-01-11    3.131125
      2010-01-12    3.035232
      2010-01-13    2.923144
      2010-01-14    2.811055
      2010-01-15    2.785824
      Freq: D

rolling_mean 함수는 약 12 ​​개 정도의 다른 함수와 함께 루 브릭 이동 창 함수 아래에 Pandas 문서에서 비공식적으로 그룹화됩니다 . Pandas의 두 번째 관련 함수 그룹은 지수 가중치 함수 (예 : 지수 이동 가중치 평균을 계산하는 ewma )라고합니다. 이 두 번째 그룹이 첫 번째 ( 이동 창 함수)에 포함되지 않는다는 사실 은 아마도 지수 가중치 변환이 고정 길이 창에 의존하지 않기 때문일 것입니다.


답변

이를 달성하는 간단한 방법은 np.convolve. 이이면의 아이디어는 이산 컨볼 루션 이 계산 되는 방식을 활용하고 이를 사용하여 롤링 평균 을 반환하는 입니다. 이것은 np.ones우리가 원하는 슬라이딩 윈도우 길이와 동일한 길이 의 시퀀스로 컨 볼빙하여 수행 할 수 있습니다 .

이를 위해 다음 함수를 정의 할 수 있습니다.

def moving_average(x, w):
    return np.convolve(x, np.ones(w), 'valid') / w

이 함수는 시퀀스의 컨볼 루션 x과 길이의 시퀀스를 취합니다 w. 선택한 mode것은 valid컨볼 루션 곱이 시퀀스가 ​​완전히 겹치는 지점에만 제공되도록하는 것입니다.


몇 가지 예 :

x = np.array([5,3,8,10,2,1,5,1,0,2])

길이 창이있는 이동 평균의 경우 다음과 2같습니다.

moving_average(x, 2)
# array([4. , 5.5, 9. , 6. , 1.5, 3. , 3. , 0.5, 1. ])

그리고 길이의 창 4:

moving_average(x, 4)
# array([6.5 , 5.75, 5.25, 4.5 , 2.25, 1.75, 2.  ])

어떻게 convolve작동합니까?

이산 컨볼 루션이 계산되는 방식에 대해 자세히 살펴 보겠습니다. 다음 함수는 np.convolve출력 값을 계산하는 방식을 복제하는 것을 목표로 합니다.

def mov_avg(x, w):
    for m in range(len(x)-(w-1)):
        yield sum(np.ones(w) * x[m:m+w]) / w 

위의 동일한 예에서 다음과 같은 결과도 얻을 수 있습니다.

list(mov_avg(x, 2))
# [4.0, 5.5, 9.0, 6.0, 1.5, 3.0, 3.0, 0.5, 1.0]

따라서 각 단계에서 수행되는 작업은 배열과 현재 사이의 내적을 취하는 것 입니다. 이 경우 np.ones(w)우리가 sum시퀀스를 직접 취하고 있다는 점을 감안할 때 곱셈 은 불필요 합니다.

Bellow는 첫 번째 출력을 계산하여 좀 더 명확하게 계산하는 방법의 예입니다. 다음과 같은 창을 원한다고 가정합니다 w=4.

[1,1,1,1]
[5,3,8,10,2,1,5,1,0,2]
= (1*5 + 1*3 + 1*8 + 1*10) / w = 6.5

그리고 다음 출력은 다음과 같이 계산됩니다.

  [1,1,1,1]
[5,3,8,10,2,1,5,1,0,2]
= (1*3 + 1*8 + 1*10 + 1*2) / w = 5.75

모든 중첩이 수행되면 시퀀스의 이동 평균을 반환합니다.


답변

이를 수행하는 다양한 방법과 몇 가지 벤치 마크가 있습니다. 가장 좋은 방법은 다른 라이브러리에서 최적화 된 코드를 사용하는 버전입니다. 이 bottleneck.move_mean방법은 아마도 가장 좋습니다. 이 scipy.convolve접근 방식은 매우 빠르고 확장 가능하며 구문 및 개념적으로 간단하지만 매우 큰 창 값에 대해서는 잘 확장되지 않습니다. numpy.cumsum당신이 순수해야 할 경우 방법은 좋은 numpy방법입니다.

참고 : 이들 중 일부 (예 bottleneck.move_mean:)는 중앙에 있지 않으며 데이터를 이동합니다.

import numpy as np
import scipy as sci
import scipy.signal as sig
import pandas as pd
import bottleneck as bn
import time as time

def rollavg_direct(a,n):
    'Direct "for" loop'
    assert n%2==1
    b = a*0.0
    for i in range(len(a)) :
        b[i]=a[max(i-n//2,0):min(i+n//2+1,len(a))].mean()
    return b

def rollavg_comprehension(a,n):
    'List comprehension'
    assert n%2==1
    r,N = int(n/2),len(a)
    return np.array([a[max(i-r,0):min(i+r+1,N)].mean() for i in range(N)])

def rollavg_convolve(a,n):
    'scipy.convolve'
    assert n%2==1
    return sci.convolve(a,np.ones(n,dtype='float')/n, 'same')[n//2:-n//2+1]

def rollavg_convolve_edges(a,n):
    'scipy.convolve, edge handling'
    assert n%2==1
    return sci.convolve(a,np.ones(n,dtype='float'), 'same')/sci.convolve(np.ones(len(a)),np.ones(n), 'same')

def rollavg_cumsum(a,n):
    'numpy.cumsum'
    assert n%2==1
    cumsum_vec = np.cumsum(np.insert(a, 0, 0))
    return (cumsum_vec[n:] - cumsum_vec[:-n]) / n

def rollavg_cumsum_edges(a,n):
    'numpy.cumsum, edge handling'
    assert n%2==1
    N = len(a)
    cumsum_vec = np.cumsum(np.insert(np.pad(a,(n-1,n-1),'constant'), 0, 0))
    d = np.hstack((np.arange(n//2+1,n),np.ones(N-n)*n,np.arange(n,n//2,-1)))
    return (cumsum_vec[n+n//2:-n//2+1] - cumsum_vec[n//2:-n-n//2]) / d

def rollavg_roll(a,n):
    'Numpy array rolling'
    assert n%2==1
    N = len(a)
    rolling_idx = np.mod((N-1)*np.arange(n)[:,None] + np.arange(N), N)
    return a[rolling_idx].mean(axis=0)[n-1:]

def rollavg_roll_edges(a,n):
    # see /programming/42101082/fast-numpy-roll
    'Numpy array rolling, edge handling'
    assert n%2==1
    a = np.pad(a,(0,n-1-n//2), 'constant')*np.ones(n)[:,None]
    m = a.shape[1]
    idx = np.mod((m-1)*np.arange(n)[:,None] + np.arange(m), m) # Rolling index
    out = a[np.arange(-n//2,n//2)[:,None], idx]
    d = np.hstack((np.arange(1,n),np.ones(m-2*n+1+n//2)*n,np.arange(n,n//2,-1)))
    return (out.sum(axis=0)/d)[n//2:]

def rollavg_pandas(a,n):
    'Pandas rolling average'
    return pd.DataFrame(a).rolling(n, center=True, min_periods=1).mean().to_numpy()

def rollavg_bottlneck(a,n):
    'bottleneck.move_mean'
    return bn.move_mean(a, window=n, min_count=1)

N = 10**6
a = np.random.rand(N)
functions = [rollavg_direct, rollavg_comprehension, rollavg_convolve,
        rollavg_convolve_edges, rollavg_cumsum, rollavg_cumsum_edges,
        rollavg_pandas, rollavg_bottlneck, rollavg_roll, rollavg_roll_edges]

print('Small window (n=3)')
%load_ext memory_profiler
for f in functions :
    print('\n'+f.__doc__+ ' : ')
    %timeit b=f(a,3)

print('\nLarge window (n=1001)')
for f in functions[0:-2] :
    print('\n'+f.__doc__+ ' : ')
    %timeit b=f(a,1001)

print('\nMemory\n')
print('Small window (n=3)')
N = 10**7
a = np.random.rand(N)
%load_ext memory_profiler
for f in functions[2:] :
    print('\n'+f.__doc__+ ' : ')
    %memit b=f(a,3)

print('\nLarge window (n=1001)')
for f in functions[2:-2] :
    print('\n'+f.__doc__+ ' : ')
    %memit b=f(a,1001)

타이밍, 작은 창 (n = 3)

Direct "for" loop :

4.14 s ± 23.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

List comprehension :
3.96 s ± 27.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

scipy.convolve :
1.07 ms ± 26.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

scipy.convolve, edge handling :
4.68 ms ± 9.69 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

numpy.cumsum :
5.31 ms ± 5.11 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

numpy.cumsum, edge handling :
8.52 ms ± 11.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Pandas rolling average :
9.85 ms ± 9.63 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

bottleneck.move_mean :
1.3 ms ± 12.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Numpy array rolling :
31.3 ms ± 91.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Numpy array rolling, edge handling :
61.1 ms ± 55.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

타이밍, 큰 창 (n = 1001)

Direct "for" loop :
4.67 s ± 34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

List comprehension :
4.46 s ± 14.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

scipy.convolve :
103 ms ± 165 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

scipy.convolve, edge handling :
272 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

numpy.cumsum :
5.19 ms ± 12.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

numpy.cumsum, edge handling :
8.7 ms ± 11.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Pandas rolling average :
9.67 ms ± 199 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

bottleneck.move_mean :
1.31 ms ± 15.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

메모리, 작은 창 (n = 3)

The memory_profiler extension is already loaded. To reload it, use:
  %reload_ext memory_profiler

scipy.convolve :
peak memory: 362.66 MiB, increment: 73.61 MiB

scipy.convolve, edge handling :
peak memory: 510.24 MiB, increment: 221.19 MiB

numpy.cumsum :
peak memory: 441.81 MiB, increment: 152.76 MiB

numpy.cumsum, edge handling :
peak memory: 518.14 MiB, increment: 228.84 MiB

Pandas rolling average :
peak memory: 449.34 MiB, increment: 160.02 MiB

bottleneck.move_mean :
peak memory: 374.17 MiB, increment: 75.54 MiB

Numpy array rolling :
peak memory: 661.29 MiB, increment: 362.65 MiB

Numpy array rolling, edge handling :
peak memory: 1111.25 MiB, increment: 812.61 MiB

메모리, 큰 창 (n = 1001)

scipy.convolve :
peak memory: 370.62 MiB, increment: 71.83 MiB

scipy.convolve, edge handling :
peak memory: 521.98 MiB, increment: 223.18 MiB

numpy.cumsum :
peak memory: 451.32 MiB, increment: 152.52 MiB

numpy.cumsum, edge handling :
peak memory: 527.51 MiB, increment: 228.71 MiB

Pandas rolling average :
peak memory: 451.25 MiB, increment: 152.50 MiB

bottleneck.move_mean :
peak memory: 374.64 MiB, increment: 75.85 MiB


답변

Pandas를 사용하는이 답변은 rolling_mean더 이상 Pandas의 일부가 아니므 로 위에서 수정되었습니다.

# the recommended syntax to import pandas
import pandas as pd
import numpy as np

# prepare some fake data:
# the date-time indices:
t = pd.date_range('1/1/2010', '12/31/2012', freq='D')

# the data:
x = np.arange(0, t.shape[0])

# combine the data & index into a Pandas 'Series' object
D = pd.Series(x, t)

이제 rolling아래의 예에서 10 일인 창 크기로 데이터 프레임 에서 함수 를 호출하면 됩니다.

d_mva10 = D.rolling(10).mean()

# d_mva is the same size as the original Series
# though obviously the first w values are NaN where w is the window size
d_mva10[:11]

2010-01-01    NaN
2010-01-02    NaN
2010-01-03    NaN
2010-01-04    NaN
2010-01-05    NaN
2010-01-06    NaN
2010-01-07    NaN
2010-01-08    NaN
2010-01-09    NaN
2010-01-10    4.5
2010-01-11    5.5
Freq: D, dtype: float64


답변

병목 현상을 사용하면 쉽게 해결할 수 있다고 생각합니다.

아래 기본 샘플을 참조하십시오.

import numpy as np
import bottleneck as bn

a = np.random.randint(4, 1000, size=(5, 7))
mm = bn.move_mean(a, window=2, min_count=1)

이것은 각 축을 따라 이동 평균을 제공합니다.

  • “mm”는 “a”의 이동 평균입니다.

  • “window”는 이동 평균을 고려할 최대 항목 수입니다.

  • “min_count”는 이동 평균을 고려할 최소 항목 수입니다 (예 : 첫 번째 요소 또는 배열에 nan 값이있는 경우).

좋은 부분은 Bottleneck이 nan 값을 처리하는 데 도움이되며 매우 효율적이라는 것입니다.


답변

경계 조건을 신중하게 처리하려는 경우 (가장자리 에서 사용 가능한 요소에서만 평균을 계산 ) 다음 함수가 트릭을 수행합니다.

import numpy as np

def running_mean(x, N):
    out = np.zeros_like(x, dtype=np.float64)
    dim_len = x.shape[0]
    for i in range(dim_len):
        if N%2 == 0:
            a, b = i - (N-1)//2, i + (N-1)//2 + 2
        else:
            a, b = i - (N-1)//2, i + (N-1)//2 + 1

        #cap indices to min and max indices
        a = max(0, a)
        b = min(dim_len, b)
        out[i] = np.mean(x[a:b])
    return out

>>> running_mean(np.array([1,2,3,4]), 2)
array([1.5, 2.5, 3.5, 4. ])

>>> running_mean(np.array([1,2,3,4]), 3)
array([1.5, 2. , 3. , 3.5])