似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
似乎没有函数可以简单地计算numpy/scipy的移动平均值,这导致了复杂的解决方案。
我的问题有两个方面:
用numpy(正确地)实现移动平均的最简单方法是什么? 既然这似乎不是小事,而且容易出错,有没有一个很好的理由不包括电池在这种情况下?
当前回答
实际上,我想要一个稍微不同于公认答案的行为。我正在为sklearn管道构建一个移动平均特征提取器,因此我要求移动平均的输出与输入具有相同的维数。我想要的是让移动平均假设级数保持不变,即[1,2,3,4,5]与窗口2的移动平均将得到[1.5,2.5,3.5,4.5,5.0]。
对于列向量(我的用例)我们得到
def moving_average_col(X, n):
z2 = np.cumsum(np.pad(X, ((n,0),(0,0)), 'constant', constant_values=0), axis=0)
z1 = np.cumsum(np.pad(X, ((0,n),(0,0)), 'constant', constant_values=X[-1]), axis=0)
return (z1-z2)[(n-1):-1]/n
对于数组
def moving_average_array(X, n):
z2 = np.cumsum(np.pad(X, (n,0), 'constant', constant_values=0))
z1 = np.cumsum(np.pad(X, (0,n), 'constant', constant_values=X[-1]))
return (z1-z2)[(n-1):-1]/n
当然,不必假设填充值为常数,但在大多数情况下这样做应该足够了。
其他回答
如果你只想要一个简单的非加权移动平均,你可以很容易地用np实现它。cumsum,可能比基于FFT的方法更快:
修正了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已经有了一些专门的功能。
如果有人需要一个简单的解决方案,这里有一个
def moving_average(a,n):
N=len(a)
return np.array([np.mean(a[i:i+n]) for i in np.arange(0,N-n+1)])
你可以通过在np.arange(0,N-n+1,step)中添加step参数来改变窗口之间的重叠
您也可以编写自己的Python C扩展。
这当然不是最简单的方法,但与使用np相比,这将使您运行得更快,内存效率更高。堆积:作为建筑块的堆积
// moving_average.c
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>
static PyObject *moving_average(PyObject *self, PyObject *args) {
PyObject *input;
int64_t window_size;
PyArg_ParseTuple(args, "Ol", &input, &window_size);
if (PyErr_Occurred()) return NULL;
if (!PyArray_Check(input) || !PyArray_ISNUMBER((PyArrayObject *)input)) {
PyErr_SetString(PyExc_TypeError, "First argument must be a numpy array with numeric dtype");
return NULL;
}
int64_t input_size = PyObject_Size(input);
double *input_data;
if (PyArray_AsCArray(&input, &input_data, (npy_intp[]){ [0] = input_size }, 1, PyArray_DescrFromType(NPY_DOUBLE)) != 0) {
PyErr_SetString(PyExc_TypeError, "Failed to simulate C array of type double");
return NULL;
}
int64_t output_size = input_size - window_size + 1;
PyObject *output = PyArray_SimpleNew(1, (npy_intp[]){ [0] = output_size }, NPY_DOUBLE);
double *output_data = PyArray_DATA((PyArrayObject *)output);
double cumsum_before = 0;
double cumsum_after = 0;
for (int i = 0; i < window_size; ++i) {
cumsum_after += input_data[i];
}
for (int i = 0; i < output_size - 1; ++i) {
output_data[i] = (cumsum_after - cumsum_before) / window_size;
cumsum_after += input_data[i + window_size];
cumsum_before += input_data[i];
}
output_data[output_size - 1] = (cumsum_after - cumsum_before) / window_size;
return output;
}
static PyMethodDef methods[] = {
{
"moving_average",
moving_average,
METH_VARARGS,
"Rolling mean of numpy array with specified window size"
},
{NULL, NULL, 0, NULL}
};
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"moving_average",
"C extension for finding the rolling mean of a numpy array",
-1,
methods
};
PyMODINIT_FUNC PyInit_moving_average(void) {
PyObject *module = PyModule_Create(&moduledef);
import_array();
return module;
}
METH_VARARGS specifies that the method only takes positional arguments. PyArg_ParseTuple allows you to parse these positional arguments. By using PyErr_SetString and returning NULL from the method, you can signal that an exception has occurred to the Python interpreter from the C extension. PyArray_AsCArray allows your method to be polymorphic when it comes to input array dtype, alignment, whether the array is C-contiguous (See "Can a numpy 1d array not be contiguous?") etc. without needing to create a copy of the array. If you instead used PyArray_DATA, you'd need to deal with this yourself. PyArray_SimpleNew allows you to create a new numpy array. This is similar to using np.empty. The array will not be initialized, and might contain non-deterministic junk which could surprise you if you forget to overwrite it.
构建C扩展
# setup.py
from setuptools import setup, Extension
import numpy
setup(
ext_modules=[
Extension(
'moving_average',
['moving_average.c'],
include_dirs=[numpy.get_include()]
)
]
)
# python setup.py build_ext --build-lib=.
基准
import numpy as np
# Our compiled C extension:
from moving_average import moving_average as moving_average_c
# Answer by Jaime using npcumsum
def moving_average_cumsum(a, n) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
# Answer by yatu using np.convolve
def moving_average_convolve(a, n):
return np.convolve(a, np.ones(n), 'valid') / n
a = np.random.rand(1_000_000)
print('window_size = 3')
%timeit moving_average_c(a, 3)
%timeit moving_average_cumsum(a, 3)
%timeit moving_average_convolve(a, 3)
print('\nwindow_size = 100')
%timeit moving_average_c(a, 100)
%timeit moving_average_cumsum(a, 100)
%timeit moving_average_convolve(a, 100)
window_size = 3
958 µs ± 4.68 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
4.52 ms ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
809 µs ± 463 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
window_size = 100
977 µs ± 937 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
6.16 ms ± 19.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
14.2 ms ± 12.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
从Numpy 1.20开始,sliding_window_view提供了一种在元素窗口中滑动/滚动的方法。然后你可以分别取平均值。
例如,对于一个4元素的窗口:
from numpy.lib.stride_tricks import sliding_window_view
# values = np.array([5, 3, 8, 10, 2, 1, 5, 1, 0, 2])
np.average(sliding_window_view(values, window_shape = 4), axis=1)
# array([6.5, 5.75, 5.25, 4.5, 2.25, 1.75, 2])
注意sliding_window_view的中间结果:
# values = np.array([5, 3, 8, 10, 2, 1, 5, 1, 0, 2])
sliding_window_view(values, window_shape = 4)
# array([[ 5, 3, 8, 10],
# [ 3, 8, 10, 2],
# [ 8, 10, 2, 1],
# [10, 2, 1, 5],
# [ 2, 1, 5, 1],
# [ 1, 5, 1, 0],
# [ 5, 1, 0, 2]])
如果你已经有一个已知大小的数组
import numpy as np
M=np.arange(12)
avg=[]
i=0
while i<len(M)-2: #for n point average len(M) - (n-1)
avg.append((M[i]+M[i+1]+M[i+2])/3) #n is denominator
i+=1
print(avg)