| 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)=
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.
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')
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.
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
Nous pouvons développer cette expression avec la fonction expand
expand_expr = expand(expr)
expand_expr
et la factoriser de nouveau sous forme factorisée avec la fonction factor
factor(expand_expr)
Nous pouvons résoudre cette expression
solve(expr)
Notez que cela équivaut à résoudre l'équation suivante pour x
[Solvers](https://docs.sympy.org/latest/modules/solvers/index.html) 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.
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
La résolution de cette équation par rapport à
solve(eq, x)
SymPy peut traiter des équations à solutions multiples
eq = Eq(expr, 1)
solve(eq, x)
La fonction solve peut également combiner plusieurs équations ensemble et résoudre un système d'équations
eq2 = Eq(x, y)
eq2
solve([eq, eq2], [x, y])
Nous pouvons aussi résoudre pour la valeur de
expr_sub = expr.subs(x, y)
expr_sub
solve(Eq(expr_sub, 1))
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
Maintenant nous simplifions cette équation en utilisant la fonction simplify
# Simplifier une expression
simplified_expr = simplify(eq)
simplified_expr
De nouveau, nous utilisons la fonction solve pour résoudre cette équation
# Résoudre l'équation
sol = solve(eq, x)
sol
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
simplify(euler)
Si cela vous intéresse, nous vous encourageons à lire le cours sur la trigonométrie et les nombres complexes.
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 :
où
Nous souhaitons calculer le point fixe de cette dynamique, c'est-à-dire la valeur de
Avec
Cela peut être facilement calculé dans SymPy
A, s, k, α, δ = symbols('A s k^* α δ')
Maintenant nous résolvons pour le point fixe
# Définir la dynamique de croissance de Solow-Swan
solow = Eq(s*A*k**α + (1-δ)*k, k)
solow
solve(solow, k)
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])
And(2*x + 5*y <= 30, x > 0)
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
Pour évaluer la somme, nous pouvons lambdify la formule.
L'expression lambdifiée peut prendre des valeurs numériques en entrée pour
sum_xy = lambdify([x, y], sum_xy)
grid = np.arange(0, 4, 1)
sum_xy(grid, grid)
Imaginez une banque avec
Elle prête
Ses dépôts sur un horizon temporel infini peuvent s'écrire
Calculons les dépôts au temps
D = symbols('D_0')
r = Symbol('r', positive=True)
Dt = Sum('(1 - r)^i * D_0', (i, 0, oo))
Dt
Nous pouvons appeler la méthode doit pour évaluer la série
Dt.doit()
Simplifier l'expression ci-dessus donne
simplify(Dt.doit())
Ceci est cohérent avec la solution dans le cours sur les séries géométriques.
Dans l'exemple suivant, nous calculons l'espérance d'une variable aléatoire discrète.
Définissons une variable aléatoire discrète
λ = symbols('lambda')
# Nous restreignons le symbole x aux entiers positifs
x = Symbol('x', integer=True, positive=True)
pmf = λ**x * exp(-λ) / factorial(x)
pmf
Nous pouvons vérifier si la somme des probabilités pour toutes les valeurs possibles égale
sum_pmf = Sum(pmf, (x, 0, oo))
sum_pmf.doit()
L'espérance de la distribution est :
fx = Sum(x*pmf, (x, 0, oo))
fx.doit()
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)
SymPy nous permet d'effectuer diverses opérations de calcul infinitésimal, telles que les limites, la différentiation et l'intégration.
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
Nous pouvons différencier toute expression SymPy en utilisant la fonction diff
# Différencier une fonction par rapport à x
df = diff(f, x)
df
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
Utilisons cette fonction pour calculer la fonction génératrice des moments de la loi exponentielle avec la fonction de densité de probabilité :
λ = 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)
Notez que nous pouvons aussi utiliser le module Stats pour calculer le moment
X = Exponential(x, λ)
moment(X, 1)
E(X**t)
En utilisant la fonction integrate, nous pouvons dériver la fonction de répartition de la loi exponentielle avec
λ_pdf = pdf.subs(λ, 1/2)
λ_pdf
integrate(λ_pdf, (x, 0, 4))
Utiliser cdf dans le module Stats donne la même solution
cdf(X, 1/2)
# Insérer une valeur pour z
λ_cdf = cdf(X, 1/2)(4)
λ_cdf
# Substituer λ
λ_cdf.subs({λ: 1/2})
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()
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()
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))
p = plot_implicit(And(2*x + 5*y <= 30, 4*x + 2*y >= 20),
(x, -1, 10), (y, -10, 10))
et des visualisations dans un espace tridimensionnel
p = plot3d(cos(2*x + y), zlabel='')
Imaginez une économie d'échange pur avec deux personnes (
Elles peuvent échanger des biens entre elles selon leurs préférences.
Supposons que les fonctions d'utilité des consommateurs sont données par
où
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
u_b
Nous nous intéressons à l'allocation Pareto-optimale des biens
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 :
# 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
# Résoudre l'équation
sol = solve(pareto, y)[0]
sol
Calculons les allocations Pareto-optimales de l'économie (courbes de contrat) avec
# Substituer α = 0.5 et β = 0.5
sol.subs({α: 0.5, β: 0.5})
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()
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?
:label: sympy_ex1
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$
:class: dropdown
Définissons d'abord la fonction
f_upper = y**x - 1
f_lower = x
f = f_upper/f_lower
f
Sympy est suffisamment intelligent pour résoudre cette limite
lim = limit(f, x, 0)
lim
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
:label: sympy_ex2
L'[estimation par maximum de vraisemblance (MLE)](https://python.quantecon.org/mle.html) 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
:class: dropdown
D'abord, nous définissons la loi binomiale
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
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
solve(Eq(log_bino_diff, 0), θ)[0]