在这一部分，我们介绍SciPy和StatsModels这两个库的使用，以及Python中的时间类型和Pandas中的时间序列。

1. SciPy

SciPy模块依赖于NumPy模块，能够提供方便快捷的N维数组操作。SciPy数组可以和NumPy的ndarray，Pandas的Series、DataFrame一起使用，并且提供包括积分、优化等许多高效的数值运算。SciPy模块不仅易于使用，而且功能强大。

在下表中，我们列出了SciPy模块主要的子模块：

1.1 线性代数

我们首先演示如何利用linalg子模块进行一些常见的线性代数运算。

首先是计算方阵的行列式：

from scipy import linalg
import numpy as np
a = np.array([[2, 3], [4, 8]])
print(linalg.det(a))

输出为：

3.999999999999999

计算方阵的逆矩阵：

a = np.array([[2, 3], [4, 8]])
b = linalg.inv(a)
print(b)

输出为：

[[ 2.   -0.75]
 [-1.    0.5 ]]

我们可以用numpy.dot来检验一下两个矩阵相乘的结果：

np.dot(a, b)

输出为一个单位阵：

array([[1., 0.],
       [0., 1.]])

我们也可以进行奇异值分解：

u, s, v = linalg.svd(a)
print(u)
print(s)
print(v)

输出为：

[[-0.37208166 -0.9282    ]
 [-0.9282      0.37208166]]
[9.63471004 0.41516558]
[[-0.46259444 -0.88657001]
 [-0.88657001  0.46259444]]

我们也可以检验一下结果，需要注意要将向量s变成一个对角线元素的矩阵：

print(u.dot(np.diag(s)).dot(v))

输出为原本的矩阵：

[[2. 3.]
 [4. 8.]]

1.2 优化和求根

我们接下来看一下optimize这个子模块的使用。

1.2.1

优化是找到最小值或等式的数值解的问题，optimize这个子模块提供了函数最小值、曲线拟合和寻找等式的根的有用算法。

我们首先定义一个函数：

def f(x):
   return x ** 2 + 10 * np.sin(x)


import matplotlib
import matplotlib.pyplot as plt
x = np.arange(-10, 10, 0.1)
plt.plot(x, f(x))
plt.show()

如下图所示：

可以看到这个函数f大概在x为-1.3时有个全局最小值,在x为3.8时有个局部最小值。接下来我们从初始点使用梯度下降法寻找最小值，例如使用BFGS算法：

from scipy import optimize
optimize.fmin_bfgs(f, 0)

打印和输出的信息如下：

Optimization terminated successfully.
         Current function value: -7.945823
         Iterations: 5
         Function evaluations: 18
         Gradient evaluations: 6
array([-1.30644012])

可以看到当初始点为0时，我们可以找到全局最优点-1.30644012，然而这个方法一个可能的问题在于，如果函数有局部最小值，会因初始点不同找到这些局部最小而不是全局最小:

optimize.fmin_bfgs(f, 5)

打印和输出的信息为：

Optimization terminated successfully.
         Current function value: 8.315586
         Iterations: 5
         Function evaluations: 18
         Gradient evaluations: 6
array([3.83746712])

例如当我们改变初始点为5时，找到的是局部最优点3.83746712。

如果我们不知道全局最小值的邻近值来选定合适的初始点，那么为了找到全局最小点，最简单的算法是蛮力算法：

# 定义要搜索的x的范围
grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
print(xmin_global)

输出为：

[-1.30641113]

当输入的搜索范围很大时，scipy.optimize.brute()会变得非常慢。在这种情况下，需要考虑更加复杂的启发式算法。

1.2.2 求根

我们可以使用scipy.optimize.fsolve来找到f(x)=0的根：

# 这里1是我们的初始猜测
root = optimize.fsolve(f, 1)
print(root)

输出为：

[0.]

我们选取的初始值同样会影响我们找到的根：

# 我们改为-3
root = optimize.fsolve(f, -3)
print(root)

输出为：

[-2.47948183]

这里我们找到了另一个根-2.47948183。

1.2.3 带约束的优化

我们最后来看如何用SciPy来求解带约束的优化问题，一般形式如下：

其中g(x)是不等式约束，h(x)是等式约束。

我们来考虑一个简单的问题，求解在x,y>=0,x+y=1的约束条件下，minimize 1-xx-yy:

from scipy.optimize import minimize
import numpy as np
e = 1e-10 # 设置一个非常接近0的值
fun = lambda x : 1 - x[0] * x[0] - x[1] * x[1]  # 目标函数
cons = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 1}, # x + y = 1
        {'type': 'ineq', 'fun': lambda x: x[0]}, # x >= 0
        {'type': 'ineq', 'fun': lambda x: x[1]}  # y >= 0
       )
x0 = np.array((1.0, 1.0)) # 设置初始值
res = minimize(fun, x0, method='SLSQP', constraints=cons)
print('最小值：',res.fun)
print('最优解：',res.x)
print('迭代终止是否成功：', res.success)
print('迭代终止原因：', res.message)

输出为：

最小值： 0.4999999999999998
最优解： [0.5 0.5]
迭代终止是否成功： True
迭代终止原因： Optimization terminated successfully.

由于这个问题的对称性，我们不难猜到最终的结果。可以看到我们使用的是minimize函数，完整的定义如下：m(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)。若是求一个函数的最大值，则改为求其相反数的最小值。上面用到的参数说明如下，更多的参数说明可以参考文档：

• fun：目标函数
• x0：初始值
• method：求解问题的算法名，默认是BFGS, L-BFGS-B, SLSQP之一
• constraints：约束条件

2. StatsModels

StatsModels模块最早起源于SciPy子模块stats中的models工具包，但后来从SciPy中移除了。到了2009年，经过修正、测试、改进并最终以全新的独立模块StatsModels诞生了。此后，StatsModels的开发团队不断添加新模型、绘图工具和统计方法，使得它最终成为了一款功能强大的统计分析工具。

作为示例，我们首先构造一下数据，并且使用OLS函数进行线性回归：

import statsmodels.api as sm
import numpy as np
x = np.linspace(0, 10, 100)
# 增加随机数
y = 3 * x + np.random.randn() + 10
# 加入一列常数项1
X = sm.add_constant(x)
# 用OLS模型进行拟合
mod = sm.OLS(y,X)
result = mod.fit()
print(result.params)

输出为：

[9.49869247 3.        ]

拟合出来的结果为y = 9.49869247 + 3x，与我们构造数据时设置的y = 10 + 3x已经非常接近。我们还可以进一步查看回归拟合的全部摘要：

print(result.summary())

输出为：

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 3.098e+32
Date:                Fri, 03 Jan 2020   Prob (F-statistic):               0.00
Time:                        15:23:48   Log-Likelihood:                 3152.7
No. Observations:                 100   AIC:                            -6301.
Df Residuals:                      98   BIC:                            -6296.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          9.4987   9.86e-16   9.63e+15      0.000       9.499       9.499
x1             3.0000    1.7e-16   1.76e+16      0.000       3.000       3.000
==============================================================================
Omnibus:                        3.156   Durbin-Watson:                   0.413
Prob(Omnibus):                  0.206   Jarque-Bera (JB):                3.146
Skew:                           0.418   Prob(JB):                        0.207
Kurtosis:                       2.762   Cond. No.                         11.7
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

我们可以使用拟合过的模型进行数据的预测：

# 预测数据
print(result.predict(X[:5]))
# 与构造的数据进行对比
print(y[:5])

输出为：

[ 9.49869247  9.80172277 10.10475308 10.40778338 10.71081368]
[ 9.49869247  9.80172277 10.10475308 10.40778338 10.71081368]

可以看到拟合的结果非常好。

StatsModels模块中另外一个重要的功能是与时间序列预测相关的ARMA模块，我们会在后面讲到线性模型的时候另行介绍。

3. 时间类型

3.1 Python中的时间类型

Python标准库包含用于日期（date）和时间（time）数据的数据类型，而且还有日历方面的功能，其中datetime以毫秒形式存储日期和时间：

from datetime import datetime
now = datetime.now()  # 直接获取现在的时间
print(now)
print(now.year)  # 年
print(now.month)  # 月
print(now.day)  # 日

输出为：

2020-01-03 12:22:06.753502
2020
1
3

timedelta表示两个datetime对象之间的时间差：

# 两个时间点之间的差：2011年1月1日零点-2008年6月24日8点15分
delta = datetime(2011, 1, 1) - datetime(2008, 6, 24, 8, 15)
print(delta)
print(delta.days)  # 一共差了多少天
print(delta.seconds)  # 一共差了多少秒

输出为：

920 days, 15:45:00
920
56700

可以给datetime对象加上（或减去）一个或多个timedelta，这样会产生一个新的datetime：

from datetime import timedelta
start = datetime(2011, 1, 7)
end = start + timedelta(12)  # 默认参数是12天
print(end)
end = start + timedelta(hours=12)  # 也可以指明别的时间单位
print(end)

输出为：

2011-01-19 00:00:00
2011-01-07 12:00:00

利用str或strftime函数（传入一个格式化字符串），datetime对象和后面要介绍的pandas的Timestamp对象可以被格式化为字符串：

stamp = datetime(2011, 1, 3)
print(str(stamp))
print(stamp.strftime('%Y-%m-%d'))  # 按照年-月-日的格式转成字符串

输出为：

2011-01-03 00:00:00
2011-01-03

datetime.strptime可以反过来用这些格式化编码将字符串转换为日期：

value = '2011-01-03'
print(datetime.strptime(value, '%Y-%m-%d'))

输出为：

2011-01-03 00:00:00

datetime.strptime是通过已知格式进行日期解析的最佳方式。但是每次都要编写格式的定义是很麻烦的事情，尤其是对于一些常见的日期格式。更简便的方式是使用dateutil这个第三方包中的parser.parse方法：

from dateutil.parser import parse
print(parse('2011-01-03'))
print(parse('Jan 31, 1997 10:45 PM'))
print(parse('6/12/2011', dayfirst=True))  # 当日出现在月前面

输出为：

2011-01-03 00:00:00
1997-01-31 22:45:00
2011-12-06 00:00:00

pandas中的to_datetime也可以解析多种不同的日期表示形式：

datestrs = ['2011-07-06 12:00:00', '2011-08-06 00:00:00']
import pandas as pd
print(pd.to_datetime(datestrs))

输出为：

DatetimeIndex(['2011-07-06 12:00:00', '2011-08-06 00:00:00'], dtype='datetime64[ns]', freq=None)

它还可以处理缺失值（None、空字符串等），产生的NaT（Not a Time）是pandas中时间戳数据的null值。

datestrs = ['', None]
print(pd.to_datetime(datestrs))

输出为：

DatetimeIndex(['NaT', 'NaT'], dtype='datetime64[ns]', freq=None)

3.2 Pandas中的时间序列

Pandas最基本的时间序列类型就是以时间戳（通常以Python字符串或datatime对象表示）为索引的Series：

import numpy as np
dates = [datetime(2011, 1, 2), datetime(2011, 1, 5),
         datetime(2011, 1, 6), datetime(2011, 1, 8),
         datetime(2011, 1, 10), datetime(2011, 1, 12)]
ts = pd.Series(np.random.randn(6), index=dates)
print(ts)

输出为：

2011-01-02   -0.412392
2011-01-05   -0.986775
2011-01-06   -0.419371
2011-01-08   -0.151056
2011-01-10   -1.332755
2011-01-12   -0.849041
dtype: float64

这些时间点实际上被放在一个DatetimeIndex中：

print(ts.index)
# pandas用NumPy的datetime64数据类型以纳秒形式存储时间戳
print(ts.index.dtype)

输出为：

DatetimeIndex(['2011-01-02', '2011-01-05', '2011-01-07', '2011-01-08',
               '2011-01-10', '2011-01-12'],
              dtype='datetime64[ns]', freq=None)
datetime64[ns]

跟其他Series一样，不同索引的时间序列之间的算术运算会自动按日期对齐：

print(ts + ts[::2])  # ts[::2] 是每隔两个数据点取一个

输出结果为：

2011-01-02   -0.824784
2011-01-05         NaN
2011-01-07   -0.838741
2011-01-08         NaN
2011-01-10   -2.665509
2011-01-12         NaN
dtype: float64

根据标签索引选取数据时，时间序列和其它的pandas.Series很像：

stamp = ts.index[2]  # 取第三个index值
print(ts[stamp])

输出为：

-0.4193706200864804

或者传入一个可以被解释为日期的字符串：

print(ts['1/10/2011'])
print(ts['20110110'])

输出为：

-1.3327546803709482
-1.3327546803709482

datetime对象也可以进行切片：

print(ts[datetime(2011, 1, 7):])  # 选取2011年1月7号及之后的时间点对应的数据

输出为：

2011-01-07   -0.419371
2011-01-08   -0.151056
2011-01-10   -1.332755
2011-01-12   -0.849041
dtype: float64

由于大部分时间序列数据都是按照时间先后排序的，因此我们也可以用不存在于该时间序列中的时间戳对其进行切片（即范围查询）：

print(ts['1/6/2011':'1/11/2011'])

输出为：

2011-01-07   -0.419371
2011-01-08   -0.151056
2011-01-10   -1.332755
dtype: float64

对于较长的时间序列，只需传入“年”或“年月”即可轻松选取数据的切片：

# 创建一个含有1000个数据点的时间序列
longer_ts = pd.Series(np.random.randn(1000), index=pd.date_range('1/1/2000', periods=1000))
print(longer_ts['2001-05'])  # 选取2001年5月的所有数据点

输出为：

2001-05-01   -0.506888
2001-05-02    0.143700
2001-05-03   -1.540290
2001-05-04   -0.624177
2001-05-05    0.559986
2001-05-06   -1.357259
2001-05-07    0.779030
2001-05-08    0.781213
2001-05-09    0.865415
2001-05-10    0.711474
2001-05-11    0.141838
2001-05-12   -0.644892
2001-05-13    2.406310
2001-05-14    1.509315
2001-05-15    0.402660
2001-05-16    0.832908
2001-05-17    0.469435
2001-05-18   -0.222889
2001-05-19    1.459801
2001-05-20    0.787919
2001-05-21   -2.410123
2001-05-22   -0.057915
2001-05-23    0.277887
2001-05-24   -0.611071
2001-05-25   -1.343405
2001-05-26    0.229819
2001-05-27   -0.504476
2001-05-28    0.011577
2001-05-29   -0.662367
2001-05-30   -0.021763
2001-05-31    0.124495
Freq: D, dtype: float64

前面提到的这些操作对DataFrame也有效。例如，对DataFrame的行进行索引：

dates = pd.date_range('1/1/2000', periods=100, freq='W-WED')  # freq取每周的周三
long_df = pd.DataFrame(np.random.randn(100, 4),
                       index=dates,
                       columns=['Colorado', 'Texas', 'New York', 'Ohio'])
print(long_df.loc['5-2001'])  # 取出2001年5月的数据

输出为：

            Colorado     Texas  New York      Ohio
2001-05-02 -2.260121 -1.379525  0.031938 -0.392022
2001-05-09 -0.474352  0.762740  1.063075  0.745917
2001-05-16 -0.438105 -0.953318 -0.081059 -0.206626
2001-05-23  0.270892 -0.395427  0.340814 -1.675902
2001-05-30 -0.418134  0.836265  0.787212 -1.698052

在某些应用场景中，可能会存在多个观测数据落在同一个时间点上的情况：

dates = pd.DatetimeIndex(['1/1/2000', '1/2/2000', '1/2/2000', '1/2/2000', '1/3/2000'])
dup_ts = pd.Series(np.arange(5), index=dates)
print(dup_ts)

输出为：

2000-01-01    0
2000-01-02    1
2000-01-02    2
2000-01-02    3
2000-01-03    4
dtype: int32

假设我们想要对具有非唯一时间戳的数据进行聚合。一个办法是使用groupby，并传入level=0：

grouped = dup_ts.groupby(level=0)  # level=0实际上是对索引进行聚合
print(grouped.mean())
print(grouped.count())

输出为：

2000-01-01    0
2000-01-02    2
2000-01-03    4
dtype: int32
2000-01-01    1
2000-01-02    3
2000-01-03    1
dtype: int64

小结

在这一部分，我们主要介绍了SciPy和StatsModels的基本使用方法，许多子模块的功能需要大家在使用的过程中进行摸索。我们还介绍了Python中的时间类型和Pandas的时间序列。