| jupytext |
|
||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| kernelspec |
|
||||||||||||||||||||||||||||||||||||||||
| translation |
|
(sympy)=
<div id="qe-notebook-header" align="right" style="text-align:right;">
<a href="https://quantecon.org/" title="quantecon.org">
<img style="width:250px;display:inline;" width="250px" src="https://assets.quantecon.org/img/qe-menubar-logo.svg" alt="QuantEcon">
</a>
</div>
与处理数值的数值库不同，SymPy 专注于直接操作数学符号和表达式。
SymPy 提供了多种功能，包括：
- 符号表达式
- 方程求解
- 化简
- 微积分
- 矩阵
- 离散数学等
这些功能使 SymPy 成为 Mathematica 等专有符号计算软件的热门开源替代品。
在本讲座中，我们将探索 SymPy 的一些功能，并演示如何使用基本的 SymPy 函数来求解经济模型。
首先导入库并初始化符号输出的打印器
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')
首先初始化一些符号以供使用
x, y, z = symbols('x y z')
符号是 SymPy 中符号计算的基本单元。
现在我们可以使用符号 x、y 和 z 来构建表达式和方程。
这里我们先构建一个简单的表达式
expr = (x+y) ** 2
expr
我们可以用 expand 函数展开这个表达式
expand_expr = expand(expr)
expand_expr
并用 factor 函数将其分解回因式分解的形式
factor(expand_expr)
我们可以求解这个表达式
solve(expr)
注意，这等价于求解以下关于 x 的方程
[Solvers](https://docs.sympy.org/latest/modules/solvers/index.html) 是一个重要的模块，包含求解不同类型方程的工具。
SymPy 中提供了多种求解器，具体取决于问题的性质。
SymPy 提供了多个函数来操作方程。
让我们用之前定义的表达式建立一个方程
eq = Eq(expr, 0)
eq
对
solve(eq, x)
SymPy 可以处理有多个解的方程
eq = Eq(expr, 1)
solve(eq, x)
solve 函数还可以将多个方程组合在一起，求解方程组
eq2 = Eq(x, y)
eq2
solve([eq, eq2], [x, y])
我们也可以通过将
expr_sub = expr.subs(x, y)
expr_sub
solve(Eq(expr_sub, 1))
下面是另一个包含符号 x 以及函数 sin、cos 和 tan 的方程示例，使用 Eq 函数构建
# 创建一个方程
eq = Eq(cos(x) / (tan(x)/sin(x)), 0)
eq
现在我们使用 simplify 函数化简这个方程
# 化简一个表达式
simplified_expr = simplify(eq)
simplified_expr
再次使用 solve 函数求解这个方程
# 求解方程
sol = solve(eq, x)
sol
SymPy 还可以处理涉及三角函数和复数的更复杂的方程。
我们使用欧拉公式来演示这一点
# 'I' 代表虚数单位 i
euler = cos(x) + I*sin(x)
euler
simplify(euler)
如果你感兴趣，我们建议你阅读关于三角函数和复数的讲座。
不动点计算在经济学和金融学中被广泛应用。
这里我们求解索洛-斯旺增长动态的不动点：
其中
我们感兴趣的是计算这一动态的不动点，即满足
在
这可以在 SymPy 中轻松计算
A, s, k, α, δ = symbols('A s k^* α δ')
现在我们求解不动点
# 定义索洛-斯旺增长动态
solow = Eq(s*A*k**α + (1-δ)*k, k)
solow
solve(solow, k)
SymPy 还允许用户定义不等式和集合运算符，并提供了丰富的操作。
reduce_inequalities([2*x + 5*y <= 30, 4*x + 2*y <= 20], [x])
And(2*x + 5*y <= 30, x > 0)
级数在经济学和统计学中被广泛使用，从资产定价到离散随机变量的期望。
我们可以使用 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
为了计算这个求和，我们可以对公式进行 lambdify 处理。
经过 lambdify 处理的表达式可以接受
sum_xy = lambdify([x, y], sum_xy)
grid = np.arange(0, 4, 1)
sum_xy(grid, grid)
假设一家银行在时间
它将
其在无限时间范围内的存款可以写成
让我们计算时间
D = symbols('D_0')
r = Symbol('r', positive=True)
Dt = Sum('(1 - r)^i * D_0', (i, 0, oo))
Dt
我们可以调用 doit 方法来计算该级数
Dt.doit()
化简上述表达式得到
simplify(Dt.doit())
这与等比数列讲座中关于部分准备金银行制度货币乘数的解一致。
在下面的例子中，我们计算离散随机变量的期望。
让我们定义一个服从泊松分布的离散随机变量
λ = symbols('lambda')
# 将符号 x 细化为正整数
x = Symbol('x', integer=True, positive=True)
pmf = λ**x * exp(-λ) / factorial(x)
pmf
我们可以验证所有可能取值的概率之和是否等于
sum_pmf = Sum(pmf, (x, 0, oo))
sum_pmf.doit()
该分布的期望为：
fx = Sum(x*pmf, (x, 0, oo))
fx.doit()
SymPy 包含一个名为 Stats 的统计子模块。
Stats 提供了内置分布和概率分布函数。
上述计算也可以使用 Stats 模块中的期望函数 E 简化为一行代码
λ = Symbol("λ", positive = True)
# 使用 sympy.stats.Poisson() 方法
X = Poisson("x", λ)
E(X)
SymPy 允许我们执行各种微积分运算，例如极限、求导和积分。
我们可以使用 limit 函数计算给定表达式的极限
# 定义一个表达式
f = x**2 / (x-1)
# 计算极限
lim = limit(f, x, 0)
lim
我们可以使用 diff 函数对任意 SymPy 表达式求导
# 对函数关于 x 求导
df = diff(f, x)
df
我们可以使用 integrate 函数计算定积分和不定积分
# 计算不定积分
indef_int = integrate(df, x)
indef_int
让我们使用这个函数来计算指数分布的矩生成函数，其概率密度函数为：
λ = Symbol('lambda', positive=True)
x = Symbol('x', positive=True)
pdf = λ * exp(-λ*x)
pdf
t = Symbol('t', positive=True)
moment_t = integrate(exp(t*x) * pdf, (x, 0, oo))
simplify(moment_t)
注意，我们也可以使用 Stats 模块来计算矩
X = Exponential(x, λ)
moment(X, 1)
E(X**t)
使用 integrate 函数，我们可以推导出
λ_pdf = pdf.subs(λ, 1/2)
λ_pdf
integrate(λ_pdf, (x, 0, 4))
使用 Stats 模块中的 cdf 函数可以得到相同的结果
cdf(X, 1/2)
# 代入 z 的值
λ_cdf = cdf(X, 1/2)(4)
λ_cdf
# 代入 λ
λ_cdf.subs({λ: 1/2})
SymPy 提供了强大的绘图功能。
首先我们使用 plot 函数绘制一个简单的函数
f = sin(2 * sin(2 * sin(2 * sin(x))))
p = plot(f, (x, -10, 10), show=False)
p.title = '简单图形'
p.show()
与 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()
它还支持绘制隐函数和不等式的可视化
p = plot_implicit(Eq((1/x + 1/y)**2, 1))
p = plot_implicit(And(2*x + 5*y <= 30, 4*x + 2*y >= 20),
(x, -1, 10), (y, -10, 10))
以及三维空间中的可视化
p = plot3d(cos(2*x + y), zlabel='')
设想一个纯交换经济，有两个人（$a$ 和
他们可以根据各自的偏好相互交换商品。
假设消费者的效用函数为
其中
首先定义符号和效用函数
# 定义符号和效用函数
x, y, α, β = symbols('x, y, α, β')
u_a = x**α * y**(1-α)
u_b = (1 - x)**β * (1 - y)**(1 - β)
u_a
u_b
我们感兴趣的是商品
注意，当在给定另一方配置的情况下，某一方的配置最优时，该点即为帕累托有效点。
用边际效用来表示：
# 当在给定另一方配置的情况下某一方的配置最优时，该点为帕累托有效点
pareto = Eq(diff(u_a, x)/diff(u_a, y),
diff(u_b, x)/diff(u_b, y))
pareto
# 求解方程
sol = solve(pareto, y)[0]
sol
让我们使用 SymPy 计算
# 代入 α = 0.5 和 β = 0.5
sol.subs({α: 0.5, β: 0.5})
我们可以使用这个结果来可视化不同参数下的更多契约曲线
# 绘制一系列 α 和 β 的取值
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()
我们邀请你调整参数，观察契约曲线如何变化，并思考以下两个问题：
- 你能想到用
numpy绘制相同图形的方法吗？ - 编写一个
numpy实现会有多困难？
:label: sympy_ex1
洛必达法则指出，对于两个函数 $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$ 时
:class: dropdown
首先定义函数
f_upper = y**x - 1
f_lower = x
f = f_upper/f_lower
f
SymPy 足够智能，可以直接求解这个极限
lim = limit(f, x, 0)
lim
我们对比洛必达法则给出的结果
lim = limit(diff(f_upper, x)/
diff(f_lower, x), x, 0)
lim
:label: sympy_ex2
[最大似然估计（MLE）](https://python.quantecon.org/mle.html) 是一种估计统计模型参数的方法。
它通常涉及最大化对数似然函数并求解一阶导数。
二项分布的概率质量函数为
$$
f(x; n, θ) = \frac{n!}{x!(n-x)!}θ^x(1-θ)^{n-x}
$$
其中 $n$ 是试验次数，$x$ 是成功次数。
假设我们观测到一系列二元结果，在 $n$ 次试验中有 $x$ 次成功。
使用 SymPy 计算 $θ$ 的最大似然估计
:class: dropdown
首先定义二项分布
n, x, θ = symbols('n x θ')
binomial_factor = (factorial(n)) / (factorial(x)*factorial(n-r))
binomial_factor
bino_dist = binomial_factor * ((θ**x)*(1-θ)**(n-x))
bino_dist
现在计算对数似然函数并求解结果
log_bino_dist = log(bino_dist)
log_bino_diff = simplify(diff(log_bino_dist, θ))
log_bino_diff
solve(Eq(log_bino_diff, 0), θ)[0]