我有一个 3D 阵列，其中包含地球表面每个网格点的海气碳通量时间序列（模型输出）。我想删除时间序列中的趋势（线性）。我遇到了这段代码：
from matplotlib import mlab
for x in xrange(40):
for y in xrange(182):
cflux_detrended[:, x, y] = mlab.detrend_linear(cflux[:, x, y])
我可以不使用 for 循环来加快速度吗？
【问题讨论】：
标签： python numpy signal-processing scipy
Scipy 有很多signal processing 工具。
使用scipy.signal.detrend() 将删除沿数据轴的线性趋势。从文档看来，整个数据集的线性趋势将从每个网格点的时间序列中减去。
import scipy.signal
cflux_detrended = scipy.signal.detrend(cflux, axis=0)
使用scipy.signal 将获得与使用原始帖子中的方法相同的结果。使用 Josef 的 detrend_separate() 函数也会返回相同的结果。
【讨论】：
这里有两个使用 numpy.linalg.lstsq 的版本。此版本使用 np.vander 创建任何多项式趋势。
警告：除示例外未测试。
我认为这样的内容将被添加到 scikits.statsmodels 中，它还没有用于去趋势的多变量版本。对于常见的趋势情况，我们可以使用 scikits.statsmodels OLS，我们也会得到所有的结果统计信息来进行估计。
# -*- coding: utf-8 -*-
"""Detrending multivariate array
Created on Fri Dec 02 15:08:42 2011
Author: Josef Perktold
http://stackoverflow.com/questions/8355197/detrending-a-time-series-of-a-multi-dimensional-array-without-the-for-loops
I should also add the multivariate version to statsmodels
"""
import numpy as np
import matplotlib.pyplot as plt
def detrend_common(y, order=1):
'''detrend multivariate series by common trend
Paramters
---------
y : ndarray
data, can be 1d or nd. if ndim is greater then 1, then observations
are along zero axis
order : int
degree of polynomial trend, 1 is linear, 0 is constant
Returns
-------
y_detrended : ndarray
detrended data in same shape as original
'''
nobs = y.shape[0]
shape = y.shape
y_ = y.ravel()
nobs_ = len(y_)
t = np.repeat(np.arange(nobs), nobs_ /float(nobs))
exog = np.vander(t, order+1)
params = np.linalg.lstsq(exog, y_)[0]
fittedvalues = np.dot(exog, params)
resid = (y_ - fittedvalues).reshape(*shape)
return resid, params
def detrend_separate(y, order=1):
'''detrend multivariate series by series specific trends
Paramters
---------
y : ndarray
data, can be 1d or nd. if ndim is greater then 1, then observations
are along zero axis
order : int
degree of polynomial trend, 1 is linear, 0 is constant
Returns
-------
y_detrended : ndarray
detrended data in same shape as original
'''
nobs = y.shape[0]
shape = y.shape
y_ = y.reshape(nobs, -1)
kvars_ = len(y_)
t = np.arange(nobs)
exog = np.vander(t, order+1)
params = np.linalg.lstsq(exog, y_)[0]
fittedvalues = np.dot(exog, params)
resid = (y_ - fittedvalues).reshape(*shape)
return resid, params
nobs = 30
sige = 0.1
y0 = 0.5 * np.random.randn(nobs,4,3)
t = np.arange(nobs)
y_observed = y0 + t[:,None,None]
for detrend_func, name in zip([detrend_common, detrend_separate],
['common', 'separate']):
y_detrended, params = detrend_func(y_observed, order=1)
print '\n\n', name
print 'params for detrending'
print params
print 'std of detrended', y_detrended.std() #should be roughly sig=0.5 (var of y0)
print 'maxabs', np.max(np.abs(y_detrended - y0))
print 'observed'
print y_observed[-1]
print 'detrended'
print y_detrended[-1]
print 'original "true"'
print y0[-1]
plt.figure()
for i in range(4):
for j in range(3):
plt.plot(y0[:,i,j], 'bo', alpha=0.75)
plt.plot(y_detrended[:,i,j], 'ro', alpha=0.75)
plt.title(name + ' detrending: blue - original, red - detrended')
plt.show()
由于 Nicholas 指出了 scipy.signal.detrend。我的去趋势分离与 scipy.signal.detrend 基本相同，但选项更少（无轴或中断）或不同（具有多项式顺序）选项。
>>> res = signal.detrend(y_observed, axis=0)
>>> (res - y0).var()
0.016931858083279336
>>> (y_detrended - y0).var()
0.01693185808327945
>>> (res - y_detrended).var()
8.402584948582852e-30
【讨论】：
我认为简单的旧列表理解是最简单的：
cflux_detrended = np.array([[mlab.detrend_linear(t) for t in kk] for kk in cflux.T])
【讨论】：