将不规则时间序列转换为相对于最近邻居的 z 分数

Question

我有一个不规则间隔索引的时间序列。我想通过减去平均值并除以每个点的标准偏差来转换数据。但是，我只想使用预定义时间距离之外的那些数据值计算均值和标准差。在下面的示例中，我使用了规则间隔的距离，但我希望它也能适应不规则的距离。

例如：

n = 20
ts = pd.Series(np.random.rand(n),
               pd.date_range('2014-05-01', periods=n, freq='T', name='Time'))

假设我想要每个点相对于该点一分钟内所有点的 zscore。

最终的结果应该类似于以下系列。

Time
2014-05-01 00:00:00    0.707107
2014-05-01 00:01:00   -0.752435
2014-05-01 00:02:00    0.866662
2014-05-01 00:03:00   -0.576136
2014-05-01 00:04:00   -0.580471
2014-05-01 00:05:00   -0.253403
2014-05-01 00:06:00   -0.076657
2014-05-01 00:07:00    1.054413
2014-05-01 00:08:00    0.095783
2014-05-01 00:09:00   -1.030982
2014-05-01 00:10:00    1.041127
2014-05-01 00:11:00   -1.028084
2014-05-01 00:12:00    0.198363
2014-05-01 00:13:00    0.851951
2014-05-01 00:14:00   -1.152701
2014-05-01 00:15:00    1.070238
2014-05-01 00:16:00   -0.395849
2014-05-01 00:17:00   -0.968585
2014-05-01 00:18:00    0.077004
2014-05-01 00:19:00    0.707107
Freq: T, dtype: float64

Answer 1

这是我一直在做的事情。请记住，这与（我怀疑您知道，否则您可能不会问这个问题）pandas rolling 功能相关但不同。对于您提供的定期间隔数据，它会很好地匹配，我们可以使用它来进行比较。

我要做的是使用np.subtract.outer 计算一个系列中所有项目与其自身的距离。

假设我们有你的时间序列ts

import pandas as pd
import numpy as np

n = 20
np.random.seed([3,1415])
data = np.random.rand(n)
tidx = pd.date_range('2014-05-01', periods=n, freq='T', name='Time')
#                                                   ^
#                                                   |
#                                            Minute Frequency
ts = pd.Series(data, tidx, name='Bliggles')

现在我可以像这样使用时间索引来计算距离

distances = pd.DataFrame(np.subtract.outer(tidx, tidx), tidx, tidx).abs()

从这里，我测试小于所需距离的距离。说这个距离叫delta

lt_delta = (distances <= delta).stack()
lt_delta = lt_delta[lt_delta]

最后，我从lt_delta 的索引中取出值，并在ts 中找到对应的值

pd.Series(ts.ix[lt_delta.index.to_series().str.get(1)].values, lt_delta.index)

我返回一个groupby 对象，所以它看起来和感觉就像调用rolling。当我将它包装在一个函数中时，它看起来像

超级功能

def groupbydelta(ts, delta):
    tidx = ts.index
    distances = pd.DataFrame(np.subtract.outer(tidx, tidx), tidx, tidx).abs()

    lt_delta = (distances <= delta).stack()
    lt_delta = lt_delta[lt_delta]
    closest = pd.Series(ts.ix[lt_delta.index.to_series().str.get(1)].values, lt_delta.index)

    return closest.groupby(level=0)

受 root 回答的启发，我编写了一个改进的 pandas/numpy 解决方案。

def groupbydelta(ts, delta):
    tidx = ts.index
    iv = pd.DataFrame({'lo': tidx - delta, 'hi': tidx + delta}, tidx)
    return pd.concat([ts.loc[r.lo:r.hi] for i, r in iv.iterrows()],
                     keys=iv.index).groupby(level=0)

让我们测试一下。我将使用delta=pd.Timedelta(1, 'm')（即一分钟）。对于我创建的时间序列，对于每个日期时间索引，我应该看到该索引、前一分钟和后一分钟。这应该等同于ts.rolling(3, center=True)，但边缘有例外。我将两者都做并比较。

gbdelta = groupbydelta(ts, pd.Timedelta(1, 'm')).mean()
rolling = ts.rolling(3, center=True).mean()

pd.concat([gbdelta, rolling], axis=1, keys=['Delta', 'Rolling']).head()

看起来不错！两者的区别在于rolling 在边缘有NaN，而gbdelta 不需要特定数量的元素，但这是设计使然。

不规则索引呢？

np.random.seed([3,1415])
n = 7200
data = np.random.rand(n)
tidx = (pd.to_datetime(['2013-02-06']) + np.random.rand(n) * pd.Timedelta(1, 'd'))
irregular_series = pd.Series(data, tidx, name='Sketch').sort_index()

并绘制irregular_series 和一些基于最近邻的过滤版本。

但你要求 zscores：

zd = (irregular_series - gbirr.mean()) / gbirr.std()

这个 z 评分有点棘手。我必须找到分组的平均值和标准差，然后将它们与原始系列一起使用。我还在想别的办法。但这已经足够流畅了。

它是什么样子的？

fig, axes = plt.subplots(1, 2, sharey=True, figsize=[10, 5])
irregular_series.plot(style='.', ax=axes[0], title='Original')
zd.plot(style='.', ax=axes[1], title='Z-Scored')

---

回答

最后，您询问了数据示例的 z 分数。为了确保我得到正确的答案...

gbd = groupbydelta(ts, pd.Timedelta(1, 'm'))

ts.sub(gbd.mean()).div(gbd.std())

Time
2014-05-01 00:00:00    0.707107
2014-05-01 00:01:00   -0.752435
2014-05-01 00:02:00    0.866662
2014-05-01 00:03:00   -0.576136
2014-05-01 00:04:00   -0.580471
2014-05-01 00:05:00   -0.253403
2014-05-01 00:06:00   -0.076657
2014-05-01 00:07:00    1.054413
2014-05-01 00:08:00    0.095783
2014-05-01 00:09:00   -1.030982
2014-05-01 00:10:00    1.041127
2014-05-01 00:11:00   -1.028084
2014-05-01 00:12:00    0.198363
2014-05-01 00:13:00    0.851951
2014-05-01 00:14:00   -1.152701
2014-05-01 00:15:00    1.070238
2014-05-01 00:16:00   -0.395849
2014-05-01 00:17:00   -0.968585
2014-05-01 00:18:00    0.077004
2014-05-01 00:19:00    0.707107
Freq: T, dtype: float64

---

时间

受 root 回答的启发，我将函数重写为基于区间。这是有道理的，它比寻找特定长度时间序列的外部差异更有效。

代码

def pirsquared(ts, delta):
    gbd = groupbydelta(ts, delta)
    return ts.sub(gbd.mean()).div(gbd.std())

cols = ['pirsquared', 'root']
ts_len = [500, 1000, 2000, 3000, 4000]
dt_len = [1, 5, 10, 20]
summary = pd.DataFrame([], pd.MultiIndex.from_product([ts_len, dt_len], names=['Points', 'Delta']), cols)
for n in ts_len:
    for d in dt_len:
        np.random.seed([3,1415])
        data = np.random.rand(n)
        tidx = (pd.to_datetime(['2013-02-06']) + np.random.rand(n) * pd.Timedelta(1, 'd'))
        ts = pd.Series(data, tidx, name='Sketch').sort_index()
        delta = pd.Timedelta(d, 'm')
        pt = timeit(lambda: pirsquared(ts, delta), number=2) / 2
        rt = timeit(lambda: root(ts, delta), number=2) / 2
        summary.loc[(n, d), cols] = pt, rt

summary.unstack().swaplevel(0, 1, 1).sort_index(1)

Answer 2

这不是pandas/numpy 解决方案，但应该提供不错的性能。本质上，要找到最近的点，您可以使用 PyPI 上的 intervaltree 包构建 Interval Tree。

intervaltree 包使用起来相当简单，并且在语法上非常类似于字典。使用此包要记住的一件事是，区间不包含上限，因此在构建树时需要填充上限。请注意，在下面的代码中，我向上限添加了额外的纳秒。

import intervaltree

def get_ts_zscore(ts, delta):
    # Get the upper and lower bounds, padding the upper bound.
    lower = ts.index - delta
    upper = ts.index + delta +  pd.Timedelta(1, 'ns')

    # Build the interval tree.
    t = intervaltree.IntervalTree().from_tuples(zip(lower, upper, ts))

    # Extract the overlaping data points for each index value.
    ts_grps = [[iv.data for iv in t[idx]]for idx in ts.index]

    # Compute the z-scores.
    ts_data = [(x - np.mean(grp))/np.std(grp, ddof=1) for x, grp in zip(ts, ts_grps)]

    return pd.Series(ts_data, ts.index)

我无法复制您的确切预期输出，可能是因为我是如何随机生成数据的？我的输出与我运行@piRSquared 的代码完全匹配，所以我很确定它是正确的。

时间

样本数据的时间安排 (n=20)：

%timeit get_ts_zscore(ts, pd.Timedelta(1, 'm'))

100 loops, best of 3: 2.89 ms per loop


%%timeit
gbd = groupbydelta(ts, pd.Timedelta(1, 'm'))
ts.sub(gbd.mean()).div(gbd.std())

100 loops, best of 3: 7.13 ms per loop

大数据的时间安排 (n=10**4)：

%timeit get_ts_zscore(ts, pd.Timedelta(1, 'm'))

1 loops, best of 3: 1.44 s per loop


%%timeit
gbd = groupbydelta(ts, pd.Timedelta(1, 'm'))
ts.sub(gbd.mean()).div(gbd.std())

1 loops, best of 3: 5.92 s per loop