import matplotlib.pyplot as plt
# 数据点
2)
np.random.seed(= 11
N = np.linspace(0, 10, N)
x = np.random.randint(10, size=(N,)).round(3)
y
'o'); plt.plot(x, y,
24 Python数值方法
24.1 插值
世界上大部分的变量都是连续的,但是我们的观测肯定是离散的。比如气温,气象站可能整点(每小时)公布一次气温,那你怎么知道1点半的气温?直觉上,你可以用1点和2点的气温的平均值,作为1点半气温的估计。这就是所谓“插值”,直接求平均(或者说线性加权)的方法属于所谓“线性插值”。
另一种情况,就是时间序列数据可能有缺失,比如多年的数据中间少了1个,也可以考虑采用插值来补全。
首先构造示例数据。假定你的数据有11个离散的观测值 \(x =0,1,2,3 ... 10\) ,绘制出图形:
24.1.1 一维插值
这里采用scipy中的interpolate模块,interpolate.interp1d()
可以提供1维插值。
参数:
- 开头2个参数是原始数据的x和y
kind
参数指定插值的方法,常用的比如线性插值linear
,二阶样条曲线quadratic
,或者直接用前值previous
等等。- 插值的方法很多,具体可见 SciPy官网。
返回值:
这个函数会根据你提供的数据以及插值的方法,返回一个新的函数f(x)
。如果你原始数据中只有1点和2点的气温,那么你调用f(1.5)
,这个函数就会告诉你它所估计1点半气温。
注意,插值的结果一定会穿过原始数据。
import numpy as np
from scipy.interpolate import interp1d
# 同样的数据,不同的插值方法可以获得不同的函数
= interp1d(x, y, kind = 'linear') # 线性
f0 = interp1d(x, y, kind = 'quadratic') # 二阶样条曲线
f1 = interp1d(x, y, kind = 'previous') # 前值
f2 = interp1d(x, y, kind = 'nearest') # 最近值 f3
虽然我们的原始数据中并没有x = 1.5的观测值,但是通过上面的这几个函数,我们可以获得x =1.5时的插值的结果。如:
# 原始数据x x
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
1.5) # 线性插值 f0(
array(7.)
1.5) # 二次样条曲线 f1(
array(7.41522187)
我们把原始数据,和插值的结果绘制一下进行对比。显而易见,线性插值可以看成把数据点直接连起来,样条曲线插值会尽可能地平滑曲线的走势。
可以看f0(1.5), f1(1.5)
在图上的位置。
注意,插值的结果一定会穿过原始数据。
= np.linspace(0,10,200) # 把0到10切得更密集
new_data_x
'o',c='k',label = '原始数据') # 原始数据点,黑色
plt.plot(x, y, = '线性插值') # 线性插值
plt.plot(new_data_x, f0(new_data_x), label = '二阶样条曲线插值') # 二阶样条曲线插值
plt.plot(new_data_x, f1(new_data_x), label ; plt.legend()
绘制原始数据和“最近值”和“前值”的插值结果。
= np.linspace(0,10,100)
new_data_x
'o',c='k' ,label = '原始数据') # 原始数据点,黑色
plt.plot(x, y, = '最近值')
plt.plot(new_data_x, f3(new_data_x),label = '前值')
plt.plot(new_data_x, f2(new_data_x),label ; plt.legend()
24.1.2 单变量样条曲线
一维插值interp1d()
有2个缺陷:
- 不能外推,即只能补全2个端点中间的值。上例中,原始数据x的区间是0-10,那么要取
f(11)
就会出错。 - 强制通过全部点,这对某些数据来说就没法插值了。
单变量样条曲线UnivariateSpline()
则可以可以解决这2点。
参数
s
:s越小,则平滑的时候考虑的数据点少,曲线看起来就比较“急转弯”,s = 0
则强制通过每一个点。默认s = None
则用所有点进行加权,曲线会尽量“平缓地转弯”。
from scipy.interpolate import UnivariateSpline
= UnivariateSpline(x,y, s = 0) # 穿过所有点
g_s0 = UnivariateSpline(x,y, s = 3)
g_s3 = UnivariateSpline(x,y) g_sn
= np.linspace(0,10.6,200) # 取值超过了原始数据
new_data_x
'o',c='k',label = '原始数据') # 原始数据点,黑色
plt.plot(x, y, = 's=0')
plt.plot(new_data_x, g_s0(new_data_x), label = 's=3')
plt.plot(new_data_x, g_s3(new_data_x), label = 's=None')
plt.plot(new_data_x, g_sn(new_data_x), label ; plt.legend()
注意,插值获得的函数可以取超过原始数据的x,曲线会按照“惯性”前进,见上图最右端x超过10的部分。
24.1.3 生成n阶段导数
注意,interp1d()
和UnivariateSpline()
,都会通过原始数据,直接返回给你一个函数f(x)
,调用这函数就可以获得曲线上的值。
其实里面有2步骤: 1. 拟合样条曲线,把数据保存在f(x)
中。 2. 求值。
有另一种做法,可以先获得样条曲线的数据tck
,我们利用这个数据,除了生成f(x)
,还可以生成其导数f'(x)
。
- 函数
splrep()
,接受原始数据x和y,返回样条曲线的信息tck
。 - 函数
splev()
,利用上一步骤生成的曲线信息tck
,求对某个x的函数值。其中参数der
为几阶导数,为0则为原函数。
那么一下两段代码是同一个意思:从原始数据x和y中生成过数据点的样条曲线,然后求x=1.5的函数值。
# 代码1
f = UnivariateSpline(x,y,s=0)
f(1.5)
# 代码2
tck = splrep(x,y,s=0)
splev(1.5, tck, der = 0)
显然两种做法曲线应该重合
from scipy.interpolate import splrep, splev
# 保存样条曲线的信息到tck
= splrep(x,y)
tck
'o',c='k',label = '原始数据')
plt.plot(x, y,
# 用splev(new_data_x, tck)求数据点
= 'splev()')
plt.plot(new_data_x, splev(new_data_x, tck), label
# 用UnivariateSpline()返回的函数求数据点
= 'US(s=0)')
plt.plot(new_data_x, g_s0(new_data_x), label ; plt.legend()
同样,用splev(der = 1)
可以获得一阶导数。
= plt.subplots(2,1,sharex=True)
fig,axes 0].plot(new_data_x, splev(new_data_x, tck, der = 0), label = 'f(x)' )
axes[1].plot(new_data_x, splev(new_data_x, tck, der = 1), label = 'f\'(x)' )
axes[0].legend()
axes[1].legend(); axes[
24.2 凸优化
这部分基本可以理解为“求极值”,包有约束和无约束的极值。
24.2.1 简单函数求极值
以前面随机生成的函数为例。
= g_s0 # 上一节生成的函数
f
= np.linspace(0,10,500)
x = f(x)
y
; plt.plot(x,y)
brute()
函数通过暴力穷举,可以求一个区间之内的函数极小值所对应的x,显然f(x)就函数的极小值。
注意这里只能求极小值,如果要求极大值,令\(g(x) = -f(x)\)即可。
import scipy.optimize as sco
# 暴力求极值
# brute() 的参数:
# 1. 要求极值的函数。
# 2. 参数区间和间隔:2到8,间隔为0.01。注意这是一个List的List,内层的每个List对应1个参数。
# 返回值:一个列表。y的极值对应的参数,每个参数一个位置,显然这里只有1个参数。
= sco.brute(f, [[2, 8, 0.01]])[0]
x1 = f(x1)
y1
= lambda x: -1 * f(x)
g = sco.brute(g, [[2, 8, 0.01]])[0]
x2 = f(x2)
y2
;
plt.plot(x,y)='blue');
plt.scatter(x1, y1, color='red'); plt.scatter(x2, y2, color
24.2.2 简单约束求极值
在约束条件下求极值:典型的就是微观经济的消费者问题:在一定的约束条件下(如收入),最大化一个函数(如效用函数)。
有两种消费品x和y,消费者具有一个CD形式的效用函数。x和y的价格为1和2,消费者的收入为m。
那么消费者的最大化问题,就可以写成:
$$ \[\begin{aligned} \max \quad u(x,y) = & x^{0.4} y^{0.6} \\ \text{s.t.} \quad x + 2y =& m \end{aligned}\]$$
= 10
m
def fun(p):
"""
要最小化的函数,即-1 * 效用函数
"""
= p
x,y return -(x**0.4 * y**0.6)
def expend(p):
"""
第一个约束,表达式应等于0
"""
= p
x,y return m - (x + 2*y)
# 参数区间
= ((0, 10), (0, 5))
bnds
# 约束
= ({'type': 'eq', 'fun': expend}) # 第1个约束是等式约束,type应为'eq'
cons
# 求极值
= sco.minimize(fun, [5, 3], bounds=bnds, constraints=cons)
result
result
fun: -3.365865430459912
jac: array([-0.33656219, -0.67320558])
message: 'Optimization terminated successfully'
nfev: 18
nit: 6
njev: 6
status: 0
success: True
x: array([4.00028954, 2.99985523])
其中result.x,就是最优化后的2个商品消费,result.fun就是乘以-1之后的最高效用。
round(result.x, 2) np.
array([4., 3.])
round(2) # 钱全部花完 expend(result.x).
0.0
如果多加一个约束,例如x是配给供应,最多消费3单位,那么问题就是:(现在不要求钱花完)
$$ \[\begin{aligned} \max \quad u(x,y) = & x^{0.4} y^{0.6} \\ \text{s.t.} \quad x + 2y \leq & m \\ x \leq & 3 \end{aligned}\]$$
def xlim(p):
"""
第一个约束,这里的表达式应大于0;
x <= 3,对应3 - x >= 0
"""
= p
x,y return 3 - x
= ({'type': 'ineq', 'fun': expend}, # 第1个约束是不等式约束,type应为'ineq'
cons 'type': 'ineq', 'fun': xlim}) # 第2个约束是不等式约束,type应为'ineq'
{
# 求极值不变
= sco.minimize(fun, [5, 3], bounds=bnds, constraints=cons)
result
result
fun: -3.290707859718669
jac: array([-0.43876103, -0.56412134])
message: 'Optimization terminated successfully'
nfev: 6
nit: 2
njev: 2
status: 0
success: True
x: array([3. , 3.5])
配给制下,最大化效用下降了。
round(2) # 钱全部花完 expend(result.x).
0.0
24.3 解方程(组)
= g_s0 # 上一节生成的函数
f
= np.linspace(0,10,500)
x = f(x)
y
;
plt.plot(x,y)0,10],[4,4]) plt.plot([
解这个函数:f(x) - 4 = 0
。
fsolve可以用于解方程:
- 第1个参数:要解的函数,必须是等于0的形式
- 第2个参数:你对解的猜测,fsolve会在你的猜测的附近进行搜索
注意:fsolve非常依赖于你的猜测,最好不要猜得太远,否则会找不到。
= lambda x: f(x) - 4 # g(x) = f(x) - 4
g
= sco.fsolve(g,[2.5, 3.5, 5, 7])
result result
array([2.35408601, 3.41291772, 5.58018654, 7.70212142])
= np.linspace(0,10,500)
x = f(x)
y
plt.plot(x,y)='red'); plt.scatter(result,f(result),color
显然,fsolve只会找到我们的猜测值附近的解。
方程组也是类似:
如2个方程解2个未知数
\[\begin{cases} 1 + y - (x-1)^2 = 0\\ 2(x+5) - y = 0 \end{cases}\]因为本案的y很容易解出来,所以比较容易绘图。但很多函数的y不能解出表达式 y = …,因此未必能绘图。
显然,让上面的方程组同时成立,必然是2个曲线的交点。
def f(x):
return (x-1)**2 - 1
def g(x):
return 2*(x+5)
= np.linspace(-10,10,100)
x
plt.plot(x,f(x)); plt.plot(x,g(x))
解为(x,y)的形式,x和y会让前面的方程组成立(结果等于0),就是我们要的解。
def equ(p):
= p # p即 (x,y)
x, y
# 两个方程,必须是等于0的形式
= 1 + y - (x-1)**2
eq1 = 2*(x+5) - y
eq2 return eq1, eq2
# 第2个参数依然是你的最佳猜测,解的形式是(x,y),显然你也要采使方程组成立的(x,y)
= sco.fsolve(equ,(-2.5,5))
result1
result1# 结果中包含了需要的(x,y)
= sco.fsolve(equ,(5,20))
result2 result2
array([ 5.74165739, 21.48331477])
result1和result2就是方程组的2个解。绘制成图形,则正好是2个交点。
plt.plot(x,f(x));
plt.plot(x,g(x))0],result1[1])
plt.scatter(result1[0],result2[1]) plt.scatter(result2[