Отклонение во временном ряду многомерного массива без циклов for

Question

У меня есть 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?

Answer 1

Сципи имеет множество инструментов для обработки сигналов инструментов. Использование scipy.signal.detrend() удалит линейный тренд вдоль оси данных. Из документации видно, что линейный тренд полного набора данных будет вычтен из временного ряда в каждой точке сетки.

import scipy.signal
cflux_detrended = scipy.signal.detrend(cflux, axis=0)

Использование scipy.signal приведет к тому же результату, что и метод в исходном посте. Использование функции Джозефа detrend_separate() также даст тот же результат.

Answer 2

Вот две версии, использующие numpy.linalg.lstsq. Эта версия использует np.vander для создания любого полиномиального тренда.

Предупреждение: не тестировалось, кроме как на примере.

Я думаю, что-то подобное будет добавлено в scikits.statsmodels, которая также не имеет многовариантной версии для трейдинга. Для случая общего тренда мы могли бы использовать OLS scikits.statsmodels, а также получить всю статистику результатов для оценки.

# -*- coding: utf-8 -*-
"""Detrending multivariate array

Created on Fri Dec 02 15:08:42 2011

Author: Josef Perktold

/6444832/otklonenie-vo-vremennom-ryadu-mnogomernogo-massiva-bez-tsiklov-for

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()

Так как Николас указал на 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

Answer 3

Я думаю, что простое понимание старого списка проще всего:

cflux_detrended = np.array([[mlab.detrend_linear(t) for t in kk] for kk in cflux.T])