Scilab (contraction de Scientific Laboratory) est un logiciel libre, développé conjointement par
l’INRIA et l’ENPC. Il est téléchargeable gratuitement à partir de
C’est un environnement de calcul numérique qui permet d’effectuer rapidement toutes les
résolutions et représentations graphiques couramment rencontrées en mathématiques
appliquées.
Scilab (qui ressemble beaucoup à Matlab) est basé sur le principe que tout calcul,
programmation ou tracé graphique peut se faire à partir de matrices rectangulaires. En Scilab,
tout est matrice : les scalaires sont des matrices, les vecteurs lignes des matrices, les vecteurs
colonnes des matrices.
1.2 Quelques commandes utiles pour commencer
1.2.1 B.A.-Ba
Dans une ligne de commande, tout ce qui suit // est ignoré, ce qui est utile pour les commentaires.
Les commandes que nous proposons sur des lignes successives sont supposées être séparées par des
retours-chariots.
A=[1,2,3;4,5,6;7,8,9],// définit une matrice 3X3
Ajouter un point virgule en fin de ligne supprime l’affichage du résultat (le calcul
est quand même effectué). Ceci évite les longs défilements à l’écran, et s’avère vite
indispensable.
x=ones(1,100);// rien n'apparaît x,// le vecteur x a bienété défini
Les résultats sont affectés par défaut à la variable ans (“answer”), qui contient donc le
résultat du dernier calcul non affecté. Toutes les variables d’une session sont globales et
conservées en mémoire. Des erreurs proviennent souvent de confusions avec des noms de
variables déjà affectés. Il faut penser à ne pas toujours utiliser les mêmes noms, ou
à libérer les variables par clear. Les variables courantes sont accessibles par who et
whos.
a=[1,2];A=[1,2;3,4];// affecte a et A 1+1,// affecte ans //who; // toutes les variables whos(),// les détails techniques clear a who;// a disparaît clear who;// a, A et ans disparaissent xbasc(),// efface le contenu de la fenêtre graphique active
La commande stacksize permet de connaître la taille de la pile utilisée par Scilab pour
stocker les variables. Si cette taille est trop faible, on peut l’ajuster grâce stacksize(n) où n est
un entier. Pour plus de détails, faire help stacksize.
L’aide en ligne est appelée par help. La commande apropos permet de retrouver les rubriques
d’aide quand on ignore le nom exact d’une fonction.
help help help apropos apropos matrix // rubriques dont le titre contient "matrix" help matrix // aide de la fonction "matrix"
1.2.2 Opérations sur les matrices
Toutes les opérations sont matricielles. Tenter une opération entre matrices de tailles non
compatibles retournera en général un message d’erreur, sauf si une des matrices est un
scalaire. Dans ce cas, l’opération (addition, multiplication, puissance) s’appliquera terme à
terme.
Les opérations terme à terme sur les matrices s’effectuent grâce à l’opérateur classique + - et
* / ^ précédé d’un point.
En résumé les différentes opérations matricielles sont : +, - addition, soustraction * , multiplication, puissance (matricielles) .* ,. multiplication terme à terme, puissance terme à terme A\b solution de A*x=b b/A solution de x*A=b ./ division terme à terme
Cette brève présentation de Scilab est inspirée des manuels ”Scilab: une introduction” de
Jean-Philippe Chancelier et ”Démarrer en Scilab” de Bernard Ycart.
xbasc(); clear n=1000; Y=rand(n,1);// Renvoie un vecteur colonne dont les entrées sont n // réalisations pseudo-aléatoires selon la loi uniforme // sur [0,1] et pseudo-indépendantes moyenne=sum(Y)/length(Y),// Moyenne arithmétique = // (Y(1)+...+Y(n))/n mediane=median(Y),// Médiane plot2d([1:n],Y),// Tracé du vecteur Y. Dans ce cas [1:n] est facultatif xset('window',1) histplot(20,Y),// Histogrammeà 20 classes
2 Loi Forte des Grands Nombres
Soit (Y1,…,Yn) un n-échantillon de variables aléatoires de loi uniforme sur [0, 1]. En regardant
l’aide sur la fonction cumsum (tapez help cumsum), calculez le vecteur Y = (Y 1,…,Y n) des
moyennes empiriques où Y i = ∑k=1iYk et tracez l’évolution de la moyenne empirique iY i à
l’aide de la fonction plot2d.
3 Théorème central limite
On considère un n-échantillon (Z1,…,Zn) où les variables aléatoires Zi sont i.i.d de même loi que
, les variables aléatoires (Ui,i ≤ p) étant i.i.d de loi uniforme sur
[0,1].
Utiliser le programme suivant pour tracer l’histogramme à nc classes de (Zi,i ≤ n). Faites
varier n, p, nc. Qu’observez-vous pour p = 1, p = 12, nc grand et nc petit? On choisira n de l’ordre
de 1000.
functiontcl(n,p,nc) xbasc(); X=rand(n,p); Z=sqrt(12/p)*(sum(X,'c')-p/2);// somme des colonnes de X, centrage et // renormalisation histplot(nc,Z) C=[-5:1/1000:5]; plot2d(C,exp(-C .^2/2)/sqrt(2*%pi),3);// densité de la loi N(0,1) endfunction
4 Cas des vecteurs gaussiens
On veut simuler des variables aléatoires indépendantes suivant la loi gaussienne centrée réduite
𝒩(0, 1).
Soient R et 𝜃 deux variables aléatoires indépendantes. On pose X = R cos(𝜃) et
Y = R sin(𝜃). Montrer que si R2 suit une loi exponentielle de paramètre et 𝜃 est
uniformément répartie sur [0, 2π], alors (X,Y ) a pour densité e−.
Montrer que si U suit la loi uniforme sur [0, 1] alors − log(U) suit la loi exponentielle
de paramètre λ.
En déduire que si (U1,U2) sont deux variables aléatoires indépendantes de loi uniforme
sur [0, 1] alors ( cos(2πU2), sin(2πU2)) est un couple de
variables aléatoires indépendantes de loi 𝒩(0, 1).
On veut maintenant simuler un vecteur gaussien centré de matrice de covariance .
Montrer que si suit la loi 𝒩 alors
suit
la loi 𝒩.
Remarque.La matrice ℒ = est en fait la décomposition de Cholesky de
M = (i.e. ℒℒT = M et ℒ triangulaire inférieure). Cette méthode permet de
simuler un vecteur gaussien dès que l’on a calculé la décomposition de Cholesky de sa
matrice de covariance. Sous Scilab la fonction chol permet de calculer la décomposition de
Cholesky d’une matrice symétrique. Cette fonction renvoie ℒT. Une version de l’algorithme
de Cholesky est proposée en annexe.
Utilisez le programme suivant pour simuler un vecteur gaussien centré et faites varier
ρ.
functiongaussian_vector(rho) if abs(rho) >1then disp('la corrélation doitêtre comprise entre -1 et 1!') disp('aborting...') return end xbasc(); n=1000; u1=rand(1,n); u2=rand(1,n); g1=ZZ;//à compléter g2=ZZ;//à compléter A=[1,0;rho,sqrt(1-rho^2)]; z=A*[g1;g2]; x=z(1,:); y=z(2,:); plot2d(x,y,-1) endfunction
5 Estimateurs à noyaux
On cherche à estimer la densité de la loi d’un échantillon. La première approche est de tracer
l’histogramme renormalisé des valeurs obtenues.
Une autre est d’estimer la densité de l’échantillon simulé par la méthode des noyaux décrite
ci-après.
Soit f la densité de probabilité à estimer. Soit (X1,...,Xn) un échantillon de la v.a. i.i.d de loi
de densité f. La mesure empirique
de
l’échantillon, où δx désigne la mesure de Dirac au point x, “converge” vers la mesure μ. Mais cette
mesure empirique n’admet pas de densité par rapport à la mesure de Lebesgue. C’est
pourquoi, on “régularise par convolution” Πn avec une suite de noyaux (Kh)h>0 qui vérifie
:
On
peut ainsi considérer la suite de noyaux Kh(x) = K(x∕h) où K peut par exemple désigner le noyau
gaussien
ou le
noyau d’Epanechnikov
On
estime alors la densité f par la fonction
La
suite n converge vers f. On admettra qu’un choix judicieux pour la suite (hn)n est de prendre
(hn)n de l’ordre de = n−1∕5.
On considère X et Y deux variables aléatoires gaussiennes centrées, réduites et
indépendantes. On définit une nouvelle variable aléatoire Z telle que Z vaut X avec
probabilité et aY + b avec probabilité où a > 0 et b ∈ ℝ. Vérifiez que la densité f de Z
s’écrit
Utilisez le programme suivant pour voir comment évolue l’estimation de la densité en fonction
de h. La vraie densité apparaît en rouge sur le graphique.
clear; // nc est le nombre de classes dans l'histogramme // n la taille de l'échantillon // h le pas de l'ordre de n^(-1/5) functionestim_noyau(n,nc,hn) xbasc(); a=1; b=3; // 2 variables alétoires uniformes U=rand(n,1); V=rand(n,1); // 2 gaussiennes centrées réduites indépendantes X=sqrt(-2*log(U)) .*cos(2*%pi*V); Y=sqrt(-2*log(U)) .*sin(2*%pi*V); // Z = X avec proba 1/3 et aY+b avec proba 2/3 epsilon=rand(n,1); Z=(a*Y+b) .*(epsilon >1/3)+X .*(epsilon <=1/3); // histogramme de la variable simulée E, nc=nombre de classes histplot(nc,Z) // estimation de la densité de la loi de E par la méthode des noyaux, // noyau Epanechnikov : C=[min(Z)-1:1/n:max(Z)+1]; for i=1:length(C)do B(i)=1/(n*hn)*3/4*sum((1-((C(i)-Z)/hn) .^2) .*(-1 <(C(i)-Z)/hn) .*((C(i)-Z)/hn <1)); end plot2d(C,B,2) //tracé de la vraie densité f=(2/(3*a)*exp(-((C-b)/a) .^2/2)+1/3*exp(-C .^2/2))/sqrt(2*%pi); plot2d(C,f,5); endfunction
6 Méthode de Monte-Carlo
La méthode de Monte-Carlo permet le calcul de valeurs approchées d’intégrales multiples en
utilisant des réalisations i.i.d. de loi uniforme. Si (Xn)n≥1 est une suite de v.a. i.i.d. de loi uniforme
sur [0, 1]m et si f : [0, 1]m→ ℝ est une fonction mesurable, alors la loi des grands nombres
appliquée à la suite de v.a.r. i.i.d. (f(Xn))n≥1 entraîne la convergence presque-sûre suivante
:
Si la
variance de f(X1) existe, le TCL assure que la convergence à lieu en 1∕. Cette vitesse est
relativement lente en petite dimension comparée aux méthodes déterministes. En revanche elle
devient très pertinente en dimension élevée. De plus cette méthode ne demande, à priori, aucune
régularité particulière sur f.
Vérifiez que les trois intégrales suivantes sont bien égales à π. Utilisez le programme suivant
pour comprendre comment approximer la valeur de π. Quelle est selon vous la plus
efficace?