Nombres aléatoires, Monte-Carlo, ...
Dans ce TD, nous allons voir comment fonctionne la méthode de
Monte-Carlo pour calculer une intégrale.
La méthode Monte-Carlo -- dont le nom fait allusion aux casinos,
temples du hasard -- utilise de grandes quantités de nombres
aléatoires. La première étape du TD est donc consacrée à la
programmation d'un générateur de nombres pseudo-aléatoires. Même s'il
existe des fonctions et dispositifs tout faits qui fournissent des
nombres aléatoires (Linux: man 3 random et man 4
random), il peut parfois être intéressant de disposer d'une
source qu'on connaît, et qui peut fournir plusieurs fois de suite
strictement la même séquence de nombres pseudo-aléatoires.
-
Programmez un générateur de nombres pseudo-aléatoires de la classe
des générateurs linéaires congruentiels:
Uj+1 = a Uj + c (mod m)
Utilisez par exemple a=4096, c=150889, m=714025. Faites de votre
générateur une fonction double aleatoire() qui renvoie
à chaque appel un nouveau nombre aléatoire entre 0 et 1. Utilisez
soit une variable globale, soit une variable locale précédée de
static pour stocker le dernier nombre généré (il est
prudent d'utiliser le type long long unsigned int, pour
être sûr de ne pas avoir d'erreur de dépassement (overflow) pendant
le calcul.
-
Tirez maintenant un grand nombre de couples (x,y) de nombres
aléatoires et comptez combien de ces points dans le plan sont à
moins d'une unité de l'origine (
x²+y²<1). Vers quoi
tend le rapport de ce nombre sur le nombre de tirages, lorsque ce
dernier tend vers l'infini ?
-
En pratique le nombre de tirages est fini. À quelle erreur
s'attend-on ? Vérifiez à l'aide d'un programme.
-
En quoi le générateur de nombres pseudo-aléatoire limite-t-il la
précision qu'on peut espérer atteindre par cette méthode
d'intégration Monte-Carlo ? Quelle sera la précision maximale du
générateur indiqué plus haut ? Testez à l'aide de votre programme.
Pour le générateur caractérisé par (a = 39373, c = 0, m =
2147483647), quelle sera la limite ?
-
Utilisez le générateur aléatoire pour générer une courbe bruitée:
f(x) = (1+bruit)*exp(-x). Générez un fichier qui comporte 50 points
dans l'intervalle
0≤x<5, avec un bruit
uniformément distribué dans un intervalle [-η,η].
-
On voudrait ajuster une fonction g(x) = exp(-ax+b) à ces points.
Écrivez la fonction Χ², somme des carrés des distances des
points de notre fichier de la fonction g(x). Calculez, en dérivant,
les conditions de moindres carrés. Réécrivez ces conditions sous
forme matricielle (rappel: vous pouvez vous aider des notes de
cours).
-
Écrivez un programme qui calcule la solution de la minimisation de
Χ². Utilisez la librairie libphynum
-
Affichez dans gnuplot les points et la fonction g(x) ajustée. Faites
l'ajustement à l'aide de gnuplot (commande
fit). Le
résultat est-il le même ?