(sympy=

23. SymPy#

23.1. Vue d’ensemble#

Contrairement aux bibliothèques numériques qui traitent des valeurs, SymPy se concentre sur la manipulation directe de symboles et d’expressions mathématiques.

SymPy fournit un large éventail de fonctionnalités comprenant

l’expression symbolique
la résolution d’équations
la simplification
le calcul infinitésimal
les matrices
les mathématiques discrètes, etc.

Ces fonctions font de SymPy une alternative open-source populaire à d’autres logiciels de calcul symbolique propriétaires tels que Mathematica.

Dans ce cours, nous explorerons certaines des fonctionnalités de SymPy et démontrerons comment utiliser les fonctions de base de SymPy pour résoudre des modèles économiques.

23.2. Pour commencer#

Importons d’abord la bibliothèque et initialisons l’imprimante pour la sortie symbolique

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

# Activer l'imprimante mathjax
init_printing(use_latex='mathjax')

23.3. Algèbre symbolique#

23.3.1. Symboles#

Nous initialisons d’abord quelques symboles avec lesquels travailler

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

Les symboles sont les unités de base du calcul symbolique dans SymPy.

23.3.2. Expressions#

Nous pouvons maintenant utiliser les symboles x, y et z pour construire des expressions et des équations.

Ici, nous construisons d’abord une expression simple

expr = (x+y) ** 2
expr

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

Nous pouvons développer cette expression avec la fonction expand

expand_expr = expand(expr)
expand_expr

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

et la factoriser de nouveau sous forme factorisée avec la fonction factor

factor(expand_expr)

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

Nous pouvons résoudre cette expression

solve(expr)

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

Notez que cela équivaut à résoudre l’équation suivante pour x

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

Note

Solvers est un module important doté d’outils pour résoudre différents types d’équations.

Une variété de solveurs sont disponibles dans SymPy selon la nature du problème.

23.3.3. Équations#

SymPy fournit plusieurs fonctions pour manipuler les équations.

Développons une équation avec l’expression que nous avons définie précédemment

eq = Eq(expr, 0)
eq

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

La résolution de cette équation par rapport à \(x\) donne le même résultat que la résolution directe de l’expression

solve(eq, x)

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

SymPy peut traiter des équations à solutions multiples

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

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

La fonction solve peut également combiner plusieurs équations ensemble et résoudre un système d’équations

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

Nous pouvons aussi résoudre pour la valeur de \(y\) en substituant simplement \(x\) par \(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]\]

Voici un autre exemple d’équation avec le symbole x et les fonctions sin, cos et tan utilisant la fonction Eq

# Créer une équation
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\]

Maintenant nous simplifions cette équation en utilisant la fonction simplify

# Simplifier une expression
simplified_expr = simplify(eq)
simplified_expr

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

De nouveau, nous utilisons la fonction solve pour résoudre cette équation

# Résoudre l'équation
sol = solve(eq, x)
sol

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

SymPy peut aussi traiter des équations plus complexes impliquant la trigonométrie et les nombres complexes.

Nous le démontrons en utilisant la formule d’Euler

# 'I' représente le nombre imaginaire i 
euler = cos(x) + I*sin(x)
euler

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

simplify(euler)

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

Si cela vous intéresse, nous vous encourageons à lire le cours sur la trigonométrie et les nombres complexes.

23.3.3.1. Exemple : calcul de point fixe#

Le calcul de point fixe est fréquemment utilisé en économie et en finance.

Ici, nous résolvons le point fixe de la dynamique de croissance de Solow-Swan :

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

où \(k_t\) est le stock de capital, \(f\) est une fonction de production, \(\delta\) est un taux de dépréciation.

Nous souhaitons calculer le point fixe de cette dynamique, c’est-à-dire la valeur de \(k\) telle que \(k_{t+1} = k_t\).

Avec \(f(k) = Ak^\alpha\), nous pouvons montrer le point fixe unique de la dynamique \(k^*\) avec papier et crayon :

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

Cela peut être facilement calculé dans SymPy

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

Maintenant nous résolvons pour le point fixe \(k^*\)

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

# Définir la dynamique de croissance de Solow-Swan
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. Inégalités et logique#

SymPy permet aussi aux utilisateurs de définir des inégalités et des opérateurs d’ensembles et fournit un large éventail d”opérations.

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. Séries#

Les séries sont largement utilisées en économie et en statistique, de l’évaluation des actifs à l’espérance de variables aléatoires discrètes.

Nous pouvons construire une simple série de sommations en utilisant la fonction Sum et les symboles 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}\]

Pour évaluer la somme, nous pouvons lambdify la formule.

L’expression lambdifiée peut prendre des valeurs numériques en entrée pour \(x\) et \(y\) et calculer le résultat

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

np.int64(36)

23.3.5.1. Exemple : dépôts bancaires#

Imaginez une banque avec \(D_0\) comme dépôt au temps \(t\).

Elle prête \((1-r)\) de ses dépôts et conserve une fraction \(r\) en réserves de trésorerie.

Ses dépôts sur un horizon temporel infini peuvent s’écrire

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

Calculons les dépôts au temps \(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}\]

Nous pouvons appeler la méthode doit pour évaluer la série

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

Simplifier l’expression ci-dessus donne

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

Ceci est cohérent avec la solution dans le cours sur les séries géométriques.

23.3.5.2. Exemple : variable aléatoire discrète#

Dans l’exemple suivant, nous calculons l’espérance d’une variable aléatoire discrète.

Définissons une variable aléatoire discrète \(X\) suivant une loi de Poisson :

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

λ = symbols('lambda')

# Nous restreignons le symbole x aux entiers positifs
x = Symbol('x', integer=True, positive=True)
pmf = λ**x * exp(-λ) / factorial(x)
pmf

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

Nous pouvons vérifier si la somme des probabilités pour toutes les valeurs possibles égale \(1\) :

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

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

\[\displaystyle 1\]

L’espérance de la distribution est :

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

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

\[\displaystyle \lambda\]

SymPy inclut un sous-module de statistique appelé Stats.

Stats offre des distributions intégrées et des fonctions sur les distributions de probabilité.

Le calcul ci-dessus peut aussi être condensé en une seule ligne en utilisant la fonction d’espérance E dans le module Stats

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

# Utilisation de la méthode sympy.stats.Poisson()
X = Poisson("x", λ)
E(X)

\[\displaystyle λ\]

23.4. Calcul symbolique#

SymPy nous permet d’effectuer diverses opérations de calcul infinitésimal, telles que les limites, la différentiation et l’intégration.

23.4.1. Limites#

Nous pouvons calculer les limites d’une expression donnée en utilisant la fonction limit

# Définir une expression
f = x**2 / (x-1)

# Calculer la limite
lim = limit(f, x, 0)
lim

\[\displaystyle 0\]

23.4.2. Dérivées#

Nous pouvons différencier toute expression SymPy en utilisant la fonction diff

# Différencier une fonction par rapport à x
df = diff(f, x)
df

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

23.4.3. Intégrales#

Nous pouvons calculer des intégrales définies et indéfinies en utilisant la fonction integrate

# Calculer l'intégrale indéfinie
indef_int = integrate(df, x)
indef_int

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

Utilisons cette fonction pour calculer la fonction génératrice des moments de la loi exponentielle avec la fonction de densité de probabilité :

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

Notez que nous pouvons aussi utiliser le module Stats pour calculer le moment

X = Exponential(x, λ)

moment(X, 1)

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

E(X**t)

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

En utilisant la fonction integrate, nous pouvons dériver la fonction de répartition de la loi exponentielle avec \(\lambda = 0.5\)

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

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

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

\[\displaystyle 0.864664716763387\]

Utiliser cdf dans le module Stats donne la même solution

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

# Insérer une valeur pour z 
λ_cdf = cdf(X, 1/2)(4)
λ_cdf

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

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

\[\displaystyle 0.864664716763387\]

23.5. Représentation graphique#

SymPy fournit une puissante fonctionnalité de représentation graphique.

D’abord, nous traçons une fonction simple en utilisant la fonction plot

f = sin(2 * sin(2 * sin(2 * sin(x))))
p = plot(f, (x, -10, 10), show=False)
p.title = 'Un graphique simple'
p.show()

_images/1179385bc9bf699824d00b14c6ca27769a83e579391e748efff89d5925bad958.png

De manière similaire à Matplotlib, SymPy fournit une interface pour personnaliser le graphique

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/2a29f1d5b2f0db1354e06555dea9138acc5d9e46050e671884b59dd2601c21d5.png

Il prend aussi en charge le tracé de fonctions implicites et la visualisation d’inégalités

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

_images/9b6f791a9506d7cc8711cc60070ba198633fbc66a00209bb2520d2c9c329ecb5.png

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

_images/deca2450156a5a78c0e0d6cf986f3e37d46c8fa1f987b0d5ad2ba312794d696b.png

et des visualisations dans un espace tridimensionnel

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

_images/4938c207047d2ac6c7cf3482968d58809282d4df0e722786afcbe26f3df59c0b.png

23.6. Application : économie d’échange à deux personnes#

Imaginez une économie d’échange pur avec deux personnes (\(a\) et \(b\)) et deux biens enregistrés en proportions (\(x\) et \(y\)).

Elles peuvent échanger des biens entre elles selon leurs préférences.

Supposons que les fonctions d’utilité des consommateurs sont données par

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

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

où \(\alpha, \beta \in (0, 1)\).

D’abord nous définissons les symboles et les fonctions d’utilité

# Définir les symboles et les fonctions d'utilité
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 - β}\]

Nous nous intéressons à l’allocation Pareto-optimale des biens \(x\) et \(y\).

Notez qu’un point est Pareto-efficace lorsque l’allocation est optimale pour une personne étant donné l’allocation pour l’autre personne.

En termes d’utilité marginale :

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

# Un point est Pareto-efficace lorsque l'allocation est optimale 
# pour une personne étant donné l'allocation pour l'autre personne

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

# Résoudre l'équation
sol = solve(pareto, y)[0]
sol

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

Calculons les allocations Pareto-optimales de l’économie (courbes de contrat) avec \(\alpha = \beta = 0.5\) en utilisant SymPy

# Substituer α = 0.5 et β = 0.5
sol.subs({α: 0.5, β: 0.5})

\[\displaystyle 1.0 x\]

Nous pouvons utiliser ce résultat pour visualiser davantage de courbes de contrat sous différents paramètres

# Tracer une plage de α et de β
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/9b32219a06ea9e39ac77ee0c2055be12aabf8e6da67018d33b095baeb0c3b162.png

Nous vous invitons à jouer avec les paramètres et à voir comment les courbes de contrat changent et à réfléchir aux deux questions suivantes :

Pouvez-vous imaginer une manière de dessiner le même graphique en utilisant numpy ?
Quelle serait la difficulté d’écrire une implémentation numpy ?

23.7. Exercices#

Exercice 23.1

La règle de L’Hôpital énonce que pour deux fonctions \(f(x)\) et \(g(x)\), si \(\lim_{x \to a} f(x) = \lim_{x \to a} g(x) = 0\) ou \(\pm \infty\), alors

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

Utilisez SymPy pour vérifier la règle de L’Hôpital pour les fonctions suivantes

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

lorsque \(x\) tend vers \(0\)

Solution

Définissons d’abord la fonction

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

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

Sympy est suffisamment intelligent pour résoudre cette limite

lim = limit(f, x, 0)
lim

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

Nous comparons le résultat suggéré par la règle de L’Hôpital

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

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

Exercice 23.2

L”estimation par maximum de vraisemblance (MLE) est une méthode pour estimer les paramètres d’un modèle statistique.

Elle implique généralement de maximiser une fonction de log-vraisemblance et de résoudre la dérivée du premier ordre.

La loi binomiale est donnée par

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

où \(n\) est le nombre d’essais et \(x\) est le nombre de succès.

Supposons que nous avons observé une série de résultats binaires avec \(x\) succès sur \(n\) essais.

Calculez le MLE de \(θ\) en utilisant SymPy

Solution

D’abord, nous définissons la loi binomiale

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

Maintenant nous calculons la fonction de log-vraisemblance et résolvons pour le résultat

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