24. SymPy#

24.1. مروری کلی#

برخلاف کتابخانه‌های عددی که با مقادیر سروکار دارند، SymPy بر دستکاری مستقیم نمادها و عبارات ریاضی تمرکز دارد.

SymPy طیف گسترده‌ای از قابلیت‌ها را فراهم می‌کند از جمله

عبارات نمادین
حل معادلات
ساده‌سازی
حساب دیفرانسیل و انتگرال
ماتریس‌ها
ریاضیات گسسته و غیره

این توابع SymPy را به یک جایگزین متن‌باز محبوب برای سایر نرم‌افزارهای محاسباتی نمادین اختصاصی مانند Mathematica تبدیل می‌کند.

در این سخنرانی، برخی از قابلیت‌های SymPy را بررسی خواهیم کرد و نحوه استفاده از توابع اولیه SymPy برای حل مدل‌های اقتصادی را نشان خواهیم داد.

24.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

# Enable the mathjax printer
init_printing(use_latex='mathjax')

24.3. جبر نمادین#

24.3.1. نمادها#

ابتدا چند نماد را برای کار با آن‌ها مقداردهی اولیه می‌کنیم

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

نمادها واحدهای اساسی برای محاسبات نمادین در SymPy هستند.

24.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 بسته به ماهیت مسئله در دسترس هستند.

24.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]\]

همچنین می‌توانیم مقدار \(y\) را با جایگزینی ساده \(x\) با \(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 آورده شده است

# Create an equation
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 ساده می‌کنیم

# Simplify an expression
simplified_expr = simplify(eq)
simplified_expr

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

دوباره از تابع solve برای حل این معادله استفاده می‌کنیم

# Solve the equation
sol = solve(eq, x)
sol

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

SymPy همچنین می‌تواند معادلات پیچیده‌تری شامل مثلثات و اعداد مختلط را مدیریت کند.

این را با استفاده از فرمول اویلر نشان می‌دهیم

# 'I' represents the imaginary number i 
euler = cos(x) + I*sin(x)
euler

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

simplify(euler)

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

اگر علاقه‌مند هستید، شما را تشویق می‌کنیم سخنرانی در مورد مثلثات و اعداد مختلط را مطالعه کنید.

24.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\) به طوری که \(k_{t+1} = k_t\).

با \(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^* \]

# Define Solow-Swan growth dynamics
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]\]

24.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\]

24.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)

24.3.5.1. مثال: سپرده‌های بانکی#

یک بانک با \(D_0\) به عنوان سپرده در زمان \(t\) را تصور کنید.

این بانک \((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}\]

این با راه‌حل در سخنرانی در مورد سری هندسی سازگار است.

24.3.5.2. مثال: متغیر تصادفی گسسته#

در مثال زیر، انتظار یک متغیر تصادفی گسسته را محاسبه می‌کنیم.

بیایید یک متغیر تصادفی گسسته \(X\) با توزیع پواسون تعریف کنیم:

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

λ = symbols('lambda')

# We refine the symbol x to positive integers
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 توزیع‌های داخلی و توابعی روی توزیع‌های احتمال ارائه می‌دهد.

محاسبه بالا همچنین می‌تواند با استفاده از تابع انتظار E در ماژول Stats به یک خط فشرده شود

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

# Using sympy.stats.Poisson() method
X = Poisson("x", λ)
E(X)

\[\displaystyle λ\]

24.4. حساب دیفرانسیل و انتگرال نمادین#

SymPy به ما اجازه می‌دهد عملیات مختلف حساب دیفرانسیل و انتگرال مانند حد، مشتق‌گیری و انتگرال‌گیری را انجام دهیم.

24.4.1. حدها#

می‌توانیم حد را برای یک عبارت داده شده با استفاده از تابع limit محاسبه کنیم

# Define an expression
f = x**2 / (x-1)

# Compute the limit
lim = limit(f, x, 0)
lim

\[\displaystyle 0\]

24.4.2. مشتقات#

می‌توانیم هر عبارت SymPy را با استفاده از تابع diff مشتق‌گیری کنیم

# Differentiate a function with respect to x
df = diff(f, x)
df

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

24.4.3. انتگرال‌ها#

می‌توانیم انتگرال‌های معین و نامعین را با استفاده از تابع integrate محاسبه کنیم

# Calculate the indefinite integral
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\]

استفاده از cdf در ماژول Stats همان راه‌حل را می‌دهد

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}\]

# Plug in a value for z 
λ_cdf = cdf(X, 1/2)(4)
λ_cdf

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

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

\[\displaystyle 0.864664716763387\]

24.5. ترسیم نمودار#

SymPy یک قابلیت ترسیم نمودار قدرتمند فراهم می‌کند.

ابتدا یک تابع ساده را با استفاده از تابع plot ترسیم می‌کنیم

f = sin(2 * sin(2 * sin(2 * sin(x))))
p = plot(f, (x, -10, 10), show=False)
p.title = 'A Simple Plot'
p.show()

_images/968dc0e94ccbd110aba1e6dd419221a5acadbdc20b0433f5ed09ab4475af7be4.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/5e23fd14eca4ef614ea864921b1d2e00f999646bdc3d7d7ff68bb8c20ce2f6c7.png

همچنین از ترسیم توابع ضمنی و تجسم نامساوی‌ها پشتیبانی می‌کند

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

_images/62ff7ecb228901f5521480663e343c0970fb76d38c4560f91d69d3da8f7a9ef3.png

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

_images/4ed6314fe9815b60bb8af352444f4fa9a3220de4632b673141d4ebbbc8e090a2.png

و تجسم‌سازی در فضای سه‌بعدی

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

_images/ab6f7f5d8ec69d83568dd5104b4ee827afe18d17d79951af743a59c056bf074c.png

24.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)\).

ابتدا نمادها و توابع مطلوبیت را تعریف می‌کنیم

# Define symbols and utility functions
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}} \]

# A point is Pareto efficient when the allocation is optimal 
# for one person given the allocation for the other person

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)}\]

# Solve the equation
sol = solve(pareto, y)[0]
sol

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

بیایید تخصیص‌های بهینه پارتو اقتصاد (منحنی‌های قرارداد) را با \(\alpha = \beta = 0.5\) با استفاده از SymPy محاسبه کنیم

# Substitute α = 0.5 and β = 0.5
sol.subs({α: 0.5, β: 0.5})

\[\displaystyle 1.0 x\]

می‌توانیم از این نتیجه برای تجسم منحنی‌های قرارداد بیشتر تحت پارامترهای مختلف استفاده کنیم

# Plot a range of αs and βs
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/2192b4c0475534fc92261e1ad681fbc0e35cdd34f7dd12a963920623b34c7f27.png

شما را دعوت می‌کنیم با پارامترها بازی کنید و ببینید منحنی‌های قرارداد چگونه تغییر می‌کنند و در مورد دو سؤال زیر فکر کنید:

آیا می‌توانید راهی برای ترسیم همان نمودار با استفاده از numpy به نظر بیاورید؟
پیاده‌سازی numpy چقدر دشوار خواهد بود؟

24.7. تمرینات#

Exercise 24.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 to Exercise 24.1

بیایید ابتدا تابع را تعریف کنیم

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 24.2

برآورد حداکثر درستنمایی (MLE) روشی برای برآورد پارامترهای یک مدل آماری است.

این روش معمولاً شامل حداکثرسازی یک تابع لگاریتم درستنمایی و حل مشتق مرتبه اول است.

توزیع دوجمله‌ای به صورت زیر داده می‌شود

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

که در آن \(n\) تعداد آزمایش‌ها و \(x\) تعداد موفقیت‌ها است.

فرض کنید ما یک سری نتایج دودویی با \(x\) موفقیت از \(n\) آزمایش مشاهده کردیم.

MLE \(θ\) را با استفاده از SymPy محاسبه کنید

Solution to Exercise 24.2

ابتدا، توزیع دوجمله‌ای را تعریف می‌کنیم

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

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

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

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

\[\displaystyle \frac{θ^{x} \left(1 - θ\right)^{n - x} n!}{x! \left(n - x\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}\]