23. SymPy#

23.1. 概述#

与处理数值的数值库不同，SymPy 专注于直接操作数学符号和表达式。

SymPy 提供了多种功能，包括：

符号表达式
方程求解
化简
微积分
矩阵
离散数学等

这些功能使 SymPy 成为 Mathematica 等专有符号计算软件的热门开源替代品。

在本讲座中，我们将探索 SymPy 的一些功能，并演示如何使用基本的 SymPy 函数来求解经济模型。

23.2. 入门#

首先导入库并初始化符号输出的打印器

from sympy import *
from sympy.plotting import plot, plot3d_parametric_line, plot3d
from sympy.solvers.inequalities import reduce_rational_inequalities
from sympy.stats import Poisson, Exponential, Binomial, density, moment, E, cdf

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl  # i18n
import matplotlib.font_manager  # i18n
FONTPATH = "_fonts/SourceHanSerifSC-SemiBold.otf"  # i18n
mpl.font_manager.fontManager.addfont(FONTPATH)  # i18n
mpl.rcParams['font.family'] = ['Source Han Serif SC']  # i18n

# 启用 mathjax 打印器
init_printing(use_latex='mathjax')

23.3. 符号代数#

23.3.1. 符号#

首先初始化一些符号以供使用

x, y, z = symbols('x y z')

符号是 SymPy 中符号计算的基本单元。

23.3.2. 表达式#

现在我们可以使用符号 x、y 和 z 来构建表达式和方程。

这里我们先构建一个简单的表达式

expr = (x+y) ** 2
expr

\[\displaystyle \left(x + y\right)^{2}\]

我们可以用 expand 函数展开这个表达式

expand_expr = expand(expr)
expand_expr

\[\displaystyle x^{2} + 2 x y + y^{2}\]

并用 factor 函数将其分解回因式分解的形式

factor(expand_expr)

\[\displaystyle \left(x + y\right)^{2}\]

我们可以求解这个表达式

solve(expr)

\[\displaystyle \left[ \left\{ x : - y\right\}\right]\]

注意，这等价于求解以下关于 x 的方程

\[ (x + y)^2 = 0 \]

Note

Solvers 是一个重要的模块，包含求解不同类型方程的工具。

SymPy 中提供了多种求解器，具体取决于问题的性质。

23.3.3. 方程#

SymPy 提供了多个函数来操作方程。

让我们用之前定义的表达式建立一个方程

eq = Eq(expr, 0)
eq

\[\displaystyle \left(x + y\right)^{2} = 0\]

对 \(x\) 求解这个方程，结果与直接求解表达式相同

solve(eq, x)

\[\displaystyle \left[ - y\right]\]

SymPy 可以处理有多个解的方程

eq = Eq(expr, 1)
solve(eq, x)

\[\displaystyle \left[ 1 - y, \ - y - 1\right]\]

solve 函数还可以将多个方程组合在一起，求解方程组

eq2 = Eq(x, y)
eq2

\[\displaystyle x = y\]

solve([eq, eq2], [x, y])

\[\displaystyle \left[ \left( - \frac{1}{2}, \ - \frac{1}{2}\right), \ \left( \frac{1}{2}, \ \frac{1}{2}\right)\right]\]

我们也可以通过将 \(x\) 替换为 \(y\) 来求 \(y\) 的值

expr_sub = expr.subs(x, y)
expr_sub

\[\displaystyle 4 y^{2}\]

solve(Eq(expr_sub, 1))

\[\displaystyle \left[ - \frac{1}{2}, \ \frac{1}{2}\right]\]

下面是另一个包含符号 x 以及函数 sin、cos 和 tan 的方程示例，使用 Eq 函数构建

# 创建一个方程
eq = Eq(cos(x) / (tan(x)/sin(x)), 0)
eq

\[\displaystyle \frac{\sin{\left(x \right)} \cos{\left(x \right)}}{\tan{\left(x \right)}} = 0\]

现在我们使用 simplify 函数化简这个方程

# 化简一个表达式
simplified_expr = simplify(eq)
simplified_expr

\[\displaystyle \cos^{2}{\left(x \right)} = 0\]

再次使用 solve 函数求解这个方程

# 求解方程
sol = solve(eq, x)
sol

\[\displaystyle \left[ - \frac{\pi}{2}, \ \frac{\pi}{2}\right]\]

SymPy 还可以处理涉及三角函数和复数的更复杂的方程。

我们使用欧拉公式来演示这一点

# 'I' 代表虚数单位 i
euler = cos(x) + I*sin(x)
euler

\[\displaystyle i \sin{\left(x \right)} + \cos{\left(x \right)}\]

simplify(euler)

\[\displaystyle e^{i x}\]

如果你感兴趣，我们建议你阅读关于三角函数和复数的讲座。

23.3.3.1. 示例：不动点计算#

不动点计算在经济学和金融学中被广泛应用。

这里我们求解索洛-斯旺增长动态的不动点：

\[ k_{t+1}=s f\left(k_t\right)+(1-\delta) k_t, \quad t=0,1, \ldots \]

其中 \(k_t\) 是资本存量，\(f\) 是生产函数，\(\delta\) 是折旧率。

我们感兴趣的是计算这一动态的不动点，即满足 \(k_{t+1} = k_t\) 的 \(k\) 值。

在 \(f(k) = Ak^\alpha\) 的条件下，我们可以用纸笔推导出唯一不动点 \(k^*\)：

\[ k^*:=\left(\frac{s A}{\delta}\right)^{1 /(1-\alpha)} \]

这可以在 SymPy 中轻松计算

A, s, k, α, δ = symbols('A s k^* α δ')

现在我们求解不动点 \(k^*\)

\[ k^* = sA(k^*)^{\alpha}+(1-\delta) k^* \]

# 定义索洛-斯旺增长动态
solow = Eq(s*A*k**α + (1-δ)*k, k)
solow

\[\displaystyle A \left(k^{*}\right)^{α} s + k^{*} \left(1 - δ\right) = k^{*}\]

solve(solow, k)

\[\displaystyle \left[ \left(\frac{A s}{δ}\right)^{- \frac{1}{α - 1}}\right]\]

23.3.4. 不等式与逻辑#

SymPy 还允许用户定义不等式和集合运算符，并提供了丰富的操作。

reduce_inequalities([2*x + 5*y <= 30, 4*x + 2*y <= 20], [x])

\[\displaystyle x \leq 5 - \frac{y}{2} \wedge x \leq 15 - \frac{5 y}{2} \wedge -\infty < x\]

And(2*x + 5*y <= 30, x > 0)

\[\displaystyle 2 x + 5 y \leq 30 \wedge x > 0\]

23.3.5. 级数#

级数在经济学和统计学中被广泛使用，从资产定价到离散随机变量的期望。

我们可以使用 Sum 函数和 Indexed 符号构建简单的求和级数

x, y, i, j = symbols("x y i j")
sum_xy = Sum(Indexed('x', i)*Indexed('y', j), 
            (i, 0, 3),
            (j, 0, 3))
sum_xy

\[\begin{split}\displaystyle \sum_{\substack{0 \leq i \leq 3\\0 \leq j \leq 3}} {x}_{i} {y}_{j}\end{split}\]

为了计算这个求和，我们可以对公式进行 lambdify 处理。

经过 lambdify 处理的表达式可以接受 \(x\) 和 \(y\) 的数值输入并计算结果

sum_xy = lambdify([x, y], sum_xy)
grid = np.arange(0, 4, 1)
sum_xy(grid, grid)

np.int64(36)

23.3.5.1. 示例：银行存款#

假设一家银行在时间 \(t\) 的存款为 \(D_0\)。

它将 \((1-r)\) 的存款贷出，并保留 \(r\) 比例的现金作为准备金。

其在无限时间范围内的存款可以写成

\[ \sum_{i=0}^\infty (1-r)^i D_0 \]

让我们计算时间 \(t\) 的存款

D = symbols('D_0')
r = Symbol('r', positive=True)
Dt = Sum('(1 - r)^i * D_0', (i, 0, oo))
Dt

\[\displaystyle \sum_{i=0}^{\infty} D_{0} \left(1 - r\right)^{i}\]

我们可以调用 doit 方法来计算该级数

Dt.doit()

\[\begin{split}\displaystyle D_{0} \left(\begin{cases} \frac{1}{r} & \text{for}\: \left|{r - 1}\right| < 1 \\\sum_{i=0}^{\infty} \left(1 - r\right)^{i} & \text{otherwise} \end{cases}\right)\end{split}\]

化简上述表达式得到

simplify(Dt.doit())

\[\begin{split}\displaystyle \begin{cases} \frac{D_{0}}{r} & \text{for}\: r > 0 \wedge r < 2 \\D_{0} \sum_{i=0}^{\infty} \left(1 - r\right)^{i} & \text{otherwise} \end{cases}\end{split}\]

这与等比数列讲座中关于部分准备金银行制度货币乘数的解一致。

23.3.5.2. 示例：离散随机变量#

在下面的例子中，我们计算离散随机变量的期望。

让我们定义一个服从泊松分布的离散随机变量 \(X\)：

\[ f(x) = \frac{\lambda^x e^{-\lambda}}{x!}, \quad x = 0, 1, 2, \ldots \]

λ = symbols('lambda')

# 将符号 x 细化为正整数
x = Symbol('x', integer=True, positive=True)
pmf = λ**x * exp(-λ) / factorial(x)
pmf

\[\displaystyle \frac{\lambda^{x} e^{- \lambda}}{x!}\]

我们可以验证所有可能取值的概率之和是否等于 \(1\)：

\[ \sum_{x=0}^{\infty} f(x) = 1 \]

sum_pmf = Sum(pmf, (x, 0, oo))
sum_pmf.doit()

\[\displaystyle 1\]

该分布的期望为：

\[ E(X) = \sum_{x=0}^{\infty} x f(x) \]

fx = Sum(x*pmf, (x, 0, oo))
fx.doit()

\[\displaystyle \lambda\]

SymPy 包含一个名为 Stats 的统计子模块。

Stats 提供了内置分布和概率分布函数。

上述计算也可以使用 Stats 模块中的期望函数 E 简化为一行代码

λ = Symbol("λ", positive = True)

# 使用 sympy.stats.Poisson() 方法
X = Poisson("x", λ)
E(X)

\[\displaystyle λ\]

23.4. 符号微积分#

SymPy 允许我们执行各种微积分运算，例如极限、求导和积分。

23.4.1. 极限#

我们可以使用 limit 函数计算给定表达式的极限

# 定义一个表达式
f = x**2 / (x-1)

# 计算极限
lim = limit(f, x, 0)
lim

\[\displaystyle 0\]

23.4.2. 导数#

我们可以使用 diff 函数对任意 SymPy 表达式求导

# 对函数关于 x 求导
df = diff(f, x)
df

\[\displaystyle - \frac{x^{2}}{\left(x - 1\right)^{2}} + \frac{2 x}{x - 1}\]

23.4.3. 积分#

我们可以使用 integrate 函数计算定积分和不定积分

# 计算不定积分
indef_int = integrate(df, x)
indef_int

\[\displaystyle x + \frac{1}{x - 1}\]

让我们使用这个函数来计算指数分布的矩生成函数，其概率密度函数为：

\[ f(x) = \lambda e^{-\lambda x}, \quad x \ge 0 \]

λ = Symbol('lambda', positive=True)
x = Symbol('x', positive=True)
pdf = λ * exp(-λ*x)
pdf

\[\displaystyle \lambda e^{- \lambda x}\]

t = Symbol('t', positive=True)
moment_t = integrate(exp(t*x) * pdf, (x, 0, oo))
simplify(moment_t)

\[\begin{split}\displaystyle \begin{cases} \frac{\lambda}{\lambda - t} & \text{for}\: \lambda > t \wedge \frac{\lambda}{t} \neq 1 \\\lambda \int\limits_{0}^{\infty} e^{x \left(- \lambda + t\right)}\, dx & \text{otherwise} \end{cases}\end{split}\]

注意，我们也可以使用 Stats 模块来计算矩

X = Exponential(x, λ)

moment(X, 1)

\[\displaystyle \frac{1}{\lambda}\]

E(X**t)

\[\displaystyle \lambda^{- t} \Gamma\left(t + 1\right)\]

使用 integrate 函数，我们可以推导出 \(\lambda = 0.5\) 时指数分布的累积密度函数

λ_pdf = pdf.subs(λ, 1/2)
λ_pdf

\[\displaystyle 0.5 e^{- 0.5 x}\]

integrate(λ_pdf, (x, 0, 4))

\[\displaystyle 0.864664716763387\]

使用 Stats 模块中的 cdf 函数可以得到相同的结果

cdf(X, 1/2)

\[\begin{split}\displaystyle \left( z \mapsto \begin{cases} 1 - e^{- z \lambda} & \text{for}\: z \geq 0 \\0 & \text{otherwise} \end{cases} \right)\end{split}\]

# 代入 z 的值
λ_cdf = cdf(X, 1/2)(4)
λ_cdf

\[\displaystyle 1 - e^{- 4 \lambda}\]

# 代入 λ
λ_cdf.subs({λ: 1/2})

\[\displaystyle 0.864664716763387\]

23.5. 绘图#

SymPy 提供了强大的绘图功能。

首先我们使用 plot 函数绘制一个简单的函数

f = sin(2 * sin(2 * sin(2 * sin(x))))
p = plot(f, (x, -10, 10), show=False)
p.title = '简单图形'
p.show()

_images/636c72526850386b38de95141bcdae87ff88c369498b5a50778832924f467fc7.png

与 Matplotlib 类似，SymPy 提供了自定义图形的接口

plot_f = plot(f, (x, -10, 10), 
              xlabel='', ylabel='', 
              legend = True, show = False)
plot_f[0].label = 'f(x)'
df = diff(f)
plot_df = plot(df, (x, -10, 10), 
            legend = True, show = False)
plot_df[0].label = 'f\'(x)'
plot_f.append(plot_df[0])
plot_f.show()

_images/167a0d236ec50b0f4d2968ebac8e4f607bbeea2336a8e09ab55e8e963b4201d3.png

它还支持绘制隐函数和不等式的可视化

p = plot_implicit(Eq((1/x + 1/y)**2, 1))

_images/0ecc4360859ccd7bc70816f969b2e5a8d308f44fac1f3553d4c2285df65bda69.png

p = plot_implicit(And(2*x + 5*y <= 30, 4*x + 2*y >= 20),
                     (x, -1, 10), (y, -10, 10))

_images/ec569101501d048de52db943967ec100bb70adf2022dcce20290138134c8e15a.png

以及三维空间中的可视化

p = plot3d(cos(2*x + y), zlabel='')

_images/a78664b1ae648ecec4e77d31f569ecee0cdc4ffed0dddf70cc72d1f8b2c22e7d.png

23.6. 应用：两人交换经济#

设想一个纯交换经济，有两个人（\(a\) 和 \(b\)）和两种商品，以比例表示（\(x\) 和 \(y\)）。

他们可以根据各自的偏好相互交换商品。

假设消费者的效用函数为

\[ u_a(x, y) = x^{\alpha} y^{1-\alpha} \]

\[ u_b(x, y) = (1 - x)^{\beta} (1 - y)^{1-\beta} \]

其中 \(\alpha, \beta \in (0, 1)\)。

首先定义符号和效用函数

# 定义符号和效用函数
x, y, α, β = symbols('x, y, α, β')
u_a = x**α * y**(1-α)
u_b = (1 - x)**β * (1 - y)**(1 - β)

u_a

\[\displaystyle x^{α} y^{1 - α}\]

u_b

\[\displaystyle \left(1 - x\right)^{β} \left(1 - y\right)^{1 - β}\]

我们感兴趣的是商品 \(x\) 和 \(y\) 的帕累托最优配置。

注意，当在给定另一方配置的情况下，某一方的配置最优时，该点即为帕累托有效点。

用边际效用来表示：

\[ \frac{\frac{\partial u_a}{\partial x}}{\frac{\partial u_a}{\partial y}} = \frac{\frac{\partial u_b}{\partial x}}{\frac{\partial u_b}{\partial y}} \]

# 当在给定另一方配置的情况下某一方的配置最优时，该点为帕累托有效点

pareto = Eq(diff(u_a, x)/diff(u_a, y), 
            diff(u_b, x)/diff(u_b, y))
pareto

\[\displaystyle \frac{y y^{1 - α} y^{α - 1} α}{x \left(1 - α\right)} = - \frac{β \left(1 - y\right) \left(1 - y\right)^{1 - β} \left(1 - y\right)^{β - 1}}{\left(1 - x\right) \left(β - 1\right)}\]

# 求解方程
sol = solve(pareto, y)[0]
sol

\[\displaystyle \frac{x β \left(α - 1\right)}{x α - x β + α β - α}\]

让我们使用 SymPy 计算 \(\alpha = \beta = 0.5\) 时经济的帕累托最优配置（契约曲线）

# 代入 α = 0.5 和 β = 0.5
sol.subs({α: 0.5, β: 0.5})

\[\displaystyle 1.0 x\]

我们可以使用这个结果来可视化不同参数下的更多契约曲线

# 绘制一系列 α 和 β 的取值
params = [{α: 0.5, β: 0.5}, 
          {α: 0.1, β: 0.9},
          {α: 0.1, β: 0.8},
          {α: 0.8, β: 0.9},
          {α: 0.4, β: 0.8}, 
          {α: 0.8, β: 0.1},
          {α: 0.9, β: 0.8},
          {α: 0.8, β: 0.4},
          {α: 0.9, β: 0.1}]

p = plot(xlabel='x', ylabel='y', show=False)

for param in params:
    p_add = plot(sol.subs(param), (x, 0, 1), 
                 show=False)
    p.append(p_add[0])
p.show()

_images/12c0c0314297869b7161a7e183c39e411619270d97516608738cd6c2f1957287.png

我们邀请你调整参数，观察契约曲线如何变化，并思考以下两个问题：

你能想到用 numpy 绘制相同图形的方法吗？
编写一个 numpy 实现会有多困难？

23.7. 练习#

Exercise 23.1

洛必达法则指出，对于两个函数 \(f(x)\) 和 \(g(x)\)，如果 \(\lim_{x \to a} f(x) = \lim_{x \to a} g(x) = 0\) 或 \(\pm \infty\)，则

\[ \lim_{x \to a} \frac{f(x)}{g(x)} = \lim_{x \to a} \frac{f'(x)}{g'(x)} \]

使用 SymPy 验证以下函数的洛必达法则

\[ f(x) = \frac{y^x - 1}{x} \]

当 \(x\) 趋近于 \(0\) 时

Solution

首先定义函数

f_upper = y**x - 1
f_lower = x
f = f_upper/f_lower
f

\[\displaystyle \frac{y^{x} - 1}{x}\]

SymPy 足够智能，可以直接求解这个极限

lim = limit(f, x, 0)
lim

\[\displaystyle \log{\left(y \right)}\]

我们对比洛必达法则给出的结果

lim = limit(diff(f_upper, x)/
            diff(f_lower, x), x, 0)
lim

\[\displaystyle \log{\left(y \right)}\]

Exercise 23.2

最大似然估计（MLE）是一种估计统计模型参数的方法。

它通常涉及最大化对数似然函数并求解一阶导数。

二项分布的概率质量函数为

\[ f(x; n, θ) = \frac{n!}{x!(n-x)!}θ^x(1-θ)^{n-x} \]

其中 \(n\) 是试验次数，\(x\) 是成功次数。

假设我们观测到一系列二元结果，在 \(n\) 次试验中有 \(x\) 次成功。

使用 SymPy 计算 \(θ\) 的最大似然估计

Solution

首先定义二项分布

n, x, θ = symbols('n x θ')

binomial_factor = (factorial(n)) / (factorial(x)*factorial(n-r))
binomial_factor

\[\displaystyle \frac{n!}{x! \left(n - r\right)!}\]

bino_dist = binomial_factor * ((θ**x)*(1-θ)**(n-x))
bino_dist

\[\displaystyle \frac{θ^{x} \left(1 - θ\right)^{n - x} n!}{x! \left(n - r\right)!}\]

现在计算对数似然函数并求解结果

log_bino_dist = log(bino_dist)

log_bino_diff = simplify(diff(log_bino_dist, θ))
log_bino_diff

\[\displaystyle θ^{- x - 1} \left(1 - θ\right)^{- n + x - 1} \left(x θ^{x} \left(1 - θ\right)^{n - x + 1} - θ^{x + 1} \left(1 - θ\right)^{n - x} \left(n - x\right)\right)\]

solve(Eq(log_bino_diff, 0), θ)[0]

\[\displaystyle \frac{x}{n}\]