Approximation de ln(x) par la méthode de Monte-Carlo

Méthode de Monte-Carlo starifiée pour le calcul de $\ln(2)$

Il s'agit de calculer une approximation de $\ln(x)$ ou de toute autre expression non polynomiale par la méthode de Monte-Carlo naïve : le principe est de « tirer » au hasard dans une cible rectangulaire et de compter le nombre de fois où la « fléchette » se plante en-dessous de la courbe représentative de la dérivée de la fonction.

Giac-XCAS
monte_carlo(primitive,xm,xM,ym,yM,nb_essais):={
    local essai,cpt,x,y,derivee;
    cpt := 0.;
    derivee := deriver(primitive);
    pour essai de 1 jusque nb_essais faire
        x := alea(xm,xM);
        y := alea(ym,yM);
        cpt := cpt + ( y <= derivee(x) ) ;
    fpour;
    retourne(cpt/nb_essais)*(xM-xm) }

On peut étudier l'écart-type de la valeur trouvée avec la valeur attendue :

Giac-XCAS
monte_carlo_sigma(primitive,xm,xM,ym,yM,nb_essais):={
   local k,liste;
   liste:=[seq( (primitive(xM) - primitive(xm)) - 
                monte_carlo(primitive,xm,xM,ym,yM,nb_essais) , k = 1..1000)];
   retourne ecart_type(liste) 
}

On obtient par exemple pour $\ln(2)$:

Giac-XCAS
monte_carlo_sigma(x->ln(x),1,2,0,1,10000)
                                  Evaluation time: 112.6
                                  0.00433957839425

On peut gagner en précision en « stratifiant » nos cibles : on découpe notre rectangle principal en petits rectangles de largeur echant

Giac-XCAS
monte_carlo_strat(primitive,xm,xM,ym,yM,nb_essais,echant):={
    h := evalf((xM-xm)/echant);
    P := NULL; app := 0;
    xh := xm;
    derivee := deriver(primitive);
    pour tranche de 1 jusque echant faire
        cpt := 0;
        pour essai de 1 jusque nb_essais faire
            x := alea(xh,xh + h);
            y := alea(ym,yM); 
            si y <= derivee(x) alors
                P := P, point([x,y],couleur=tranche);
                cpt := cpt +1;
            fsi;
        fpour;
        app := app + (cpt/nb_essais)*h;
        xh := xh + h;
    fpour; 
    afficher(primitive(xM)-primitive(xm) + " vaut environ " + app);
    afficher("precision :" + (primitive(xM)-primitive(xm) - app));
    retourne(P) 
}

Par exemple, on obtient une valeur approchée de $\ln(2)$ à $2\times 10^{-4}$ près :

Giac-XCAS
monte_carlo_strat(x -> ln(x),1,2,0,1,10000,10)

courtesy of webmatter.de