用Python（NumPy或Pandas）进行时间序列的线性回归

Question

我是一个新手，对Python和编程一般性知识都不熟悉，请原谅我可能会犯的简单错误或显而易见的错误。

我的目标很简单，就是想将线性趋势(1-d多项式)拟合到一堆时间序列上，以查看斜率是正还是负。目前，我只想让它在一个时间序列上工作。

问题在于：似乎pandas和numpy都无法进行日期回归分析。我的日期并不固定（通常每个月1天但不是同一天），因此不能使用Linear Regression from Time Series Pandas中提出的建议。

我的时间序列CSV文件看起来像：

StationName,    year,   month,  day,    depth,  NO3-N,  PO4-P,  TotP,   TotN,
Kvarnbacken (Savaran),  2003,   2,  25, 0.5,    46, 9,  14, 451
Kvarnbacken (Savaran),  2003,   3,  18, 0.5,    64, 15, 17, 310
Kvarnbacken (Savaran),  2003,   3,  31, 0.5,    76, 7,  19, 566

到目前为止，我所拥有的是：

import datetime as dt
from scipy import stats
import numpy as np

# read in station csv file
data = pd.read_csv('Kvarnbacken (Savaran)_2003.csv')
data.head()
# set up dates to something python can recognize
data['date'] = pd.to_datetime(data.year*10000+data.month *
                          100+data.day, format='%Y%m%d')

我尝试了

slope, intercept, r_value, p_value, std_err = stats.linregress(data.date,
                                                               data.TotP)

and got the error TypeError: ufunc add cannot use operands with types dtype('

I also tried

coefP = np.polyfit(data.date, data.TotP, 1)
polyP = np.poly1d(coefP)
ys = polyP(data.date)
print 'For P: coef, poly'
print coefP
print polyP

and got the same error.

I am guessing the easiest way around this is to do something where I just count the days since the first measurement I have and then just do a regression with days_since to the total phosphorous concentration (totP) but I am not sure of the easiest way to do that or if there was another trick.

- I. F.

---

我认为你所做的似乎是错误的，因为两天之间的间隔不是恒定的（例如：在12月31日和1月1日之间）。你可以将第一天称为第0天。然后从每个日期中减去这个第一天，并将其转换为天数。 - Hugo

雨果，你绝对是正确的。我没有认真考虑我的快速修复方案。 - I. F.

注释＃1不是一个好的方法。如果你想做类似的事情，你应该把它转换成一些公共单位。更像365年+30月+天。尽管即使这样也不是理想的，因为年份和月份没有恒定的天数。请参见答案以了解更好的方法。 - JohnE

Answer 1

你可以按以下方式将日期时间转换为天数。

data['days_since'] = (data.date - pd.to_datetime('2003-02-25') ).astype('timedelta64[D]')

        date  days_since
0 2003-02-25           0
1 2003-03-18          21
2 2003-03-31          34

现在你应该能够像之前那样进行回归。

slope, intercept, r_value, p_value, std_err = stats.linregress(data.days_since, 
                                                               data.TotP)
slope, intercept
(0.1466591166477916, 13.977916194790488)

你可能还想考虑其他回归选项，例如statsmodels软件包，特别是如果你将经常进行这种操作。（请注意，与linregress相比，x和y被颠倒了）

import statsmodels.formula.api as smf

smf.ols( 'TotP ~ days_since', data=data ).fit().params

Intercept     13.977916
days_since     0.146659

顺便说一下，这只是statsmodels输出的一小部分（使用summary()代替params可以获得额外的输出）。

- JohnE

---

谢谢。我最终做了类似的事情，但是在如何将'day_since'转换为整数方面遇到了困难，astype是个不错的技巧。 - I. F.