SciPy和StatsModels的使用

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

1. SciPy

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

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

子模块名称 功能
cluster 聚类算法
constants 物理数学常数
fftpack 快速傅里叶变换
integrate 积分和常微分方程求解
interpolate 插值处理
io 输入输出
linalg 线性代数
ndimage 图像处理模块
odr 正交距离回归
optimize 优化和求根
signal 信号处理
sparse 稀疏矩阵
spatial 空间数据结构和算法
special 特殊函数模块
stats 统计分布和函数

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                            

剩余60%内容,订阅专栏后可继续查看/也可单篇购买

<p> 在这个专刊中,我们会学习Python在金融中的应用: ·掌握Python进行数据分析和处理的方法; ·了解Python在量化投资中的应用; ·转行到金融数据分析和量化交易领域的基础等内容。 这个专刊适合: ·想要学习Python编程的在校学生; ·想要转行到金融领域的在职人士; ·想要利用业余时间进行量化投资的在职人士等。 本专刊购买后即可解锁所有章节,故不可以退换哦~ </p>

全部评论

相关推荐

点赞 评论 收藏
分享
03-30 19:30
石家庄学院 Java
野蛮的柯基在游泳:都能入股了,还得是Java
点赞 评论 收藏
分享
评论
点赞
收藏
分享

创作者周榜

更多
牛客网
牛客企业服务