Simuler la loi exponentielle sans fonction exp ou log

Comment simuler une variable aléatoire obéissant à une loi exponentielle (du type durée de vie d'une particule radioactive), sans jamais utiliser dans le code la fonction exp ni la fonction log ?
Règle du jeu de ce fil : expliquer clairement ce qu'on fait, publier les codes.

Bonjour,

Comme à son habitude GBZM radote, il a déjà posé cette question plusieurs dizaines de fois...

Un peu de contexte, les plus jeunes d'entre vous (ceux qui ont moins de 50 ans) ne voit pas l'intéret de la question de GBZM, mais à l'époque (année 70-80) c'était gourmand en CPU de calculer des exponentielles et autres fonctions complexes (ce qui n'est pas le cas maintenant, à cause de la puissance de calcul plus importantes).

Pour éviter de se faire embobiner par GBZM qui voient du classique partout, il n'y a pas d'approximation connu* de loi expo par des binomiales (je ne parle bien sûr pas de loi de Poisson, où cela est effectivement classique).

@Dlzlogic : il est évident qu'il essait de te tirer les vers du nez, même aprés avoir dit à plusieurs reprises, que ta façon d'approximer la loi expo, ne l'intéresse pas. Comme à son habitude quand tu lui auras expliquer comment tu procédes, il risque de minimiser ton travail...

* : https://fr.wikipedia.org/wiki/Loi_binomiale

Bonne journée.

Bonjour,
Petite précision indispensable pour Hunx3 :
Dans tous les livres de maths avec un chapitre probabilité, il y a un paragraphe "loi exponentielle". Vous appliquez ce chapitre.
Ce n'est pas du tout la question.
Je la repose. "Je" ne connais pas la loi exponentielle, ni la loi géométrique, ni la loi binomiale, mais je fais de la recherche concernant le stock de led à prévoir et je cherche un modèle qui me permettra de répondre à la question.
Comme j'ai pas le temps d'attendre pour établir mon modèle, j'utilise mon ordinateur et mes capacités intellectuelles pour résoudre le problème.

PS Le premier message du fil est amusant, puisque j'ai fait cette expérience mais je refuse, pour l'instant, de donner des détails de ma méthode.
Traduction : Humx3 utilise les techniques utilisées par certains étudiants du type "je sais le faire mais j'aimerais savoir ce que vous en pensez", c'est souvent dans le chapitre "défis et énigmes".

Bien, voici comment je procède. C'est tout à fait élémentaire.
La première procédure produit une liste de temps de vie. Pour cela on choisit une probabilité p de disparaître par unité de temps (on prendra p petit). Le paramètre n est le nombre de temps de vie qu'on produit.

Code:

import numpy as np
import random as rd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def liste_temps_vie(p,n) :
    L=[]
    for _ in range(n) :
        t=0 ; vie=0
        while vie==0 :
            t+=1
            if rd.random() < p :
                vie=1
        L.append(t)
    return L

On produit une liste de 10 000 temps de vie pour p=1/1000

Code:

%time TV=liste_temps_vie(1/1000,10000)

CPU times: user 2.7 s, sys: 0 ns, total: 2.7 s
Wall time: 2.73 s

Que faire de cette liste ? On peut déjà regarder sa moyenne et sa médiane.

Code:

print("moyenne : {:.1f} ; médiane : {}".format(np.mean(TV),np.median(TV)))

moyenne : 996.4 ; médiane : 689.0

L'espérance du temps de vie est 1000, puisque p=1/1000. La moyenne tend vers l'espérance quand la taille de l'échantillon tend vers l'infini, avec un échantillon de taille 10 000 on n'est pas mauvais. On se souvient que ln(2) = 0.693...  . Le rapport médiane/moyenne en est vraiment très proche.

La suite à venir.

Pour éviter à GBZM de se répeter et répeter, et répéter, et....

Même s'il semble beaucoup aimer cela :

HumHumHum a écrit:

Code:

def distr_exp_approx(n) :
    nb=0 ; succ=0
    while succ==0 :
        nb+=1
        if rd.random() < 1/n :
            succ=1
    return nb/n

C'est très simple : cette fonction prend en entrée l'entier n. Elle fait des épreuves de Bernoulli de probabilité de succès 1/n jusqu'à ce qu'un succès (= désintégration) soit obtenu. Elle retourne le nombre d'épreuves, divisé par n : c'est la durée de vie.
On peut faire une variante, spécifiquement pour la durée de vie d'un noyau de carbone 14. Chaque épreuve de Bernoulli a probabilité de succès (= désintégration) de 1.21 * 10^-3, et on retourne le nombre d'épreuves pour obtenir un succès multiplié par 10

Code:

def vie_C14() :
    nb=0 ; succ=0
    while succ==0 :
        nb+=1
        if rd.random() < 0.00121 :
            succ=1
    return 10*nb

On peut utiliser cette fonction pour faire afficher un bel histogramme :

Code:

fig, ax = plt.subplots(1, 1)
r = np.array([vie_C14() for _ in range(10000)])
ax.hist(r,  bins='auto', histtype='stepfilled', alpha=0.5,\
        range=(0,20000))
plt.title("Durée de vie du C14 pour 10 000 noyaux au départ")
plt.show()

HumHumHum a écrit:
Du genre de ce que j'avais fait avec la distribution du nombre de noyaux de carbone-14 survivants parmi 1000 au bout d'un nombre donné d'années, c'est cela ? Ça me semble plus intéressant que cette histoire bien connue de tirage à pile ou face :

Pour le nombre de piles sur 1000 tirages à pile ou face, ça ne sera pas différent de l'histogramme du milieu.

Merci Dattier pour le grand intérêt que vous semblez porter à ce sujet.
Si vous avez une idée originale, je vous en prie, contribuez plutôt que de polluer le fil. La deuxième citation que vous faites est hors-sujet, c'est une illustration du TCL.

Aller, je te laisse tourner en rond.

Reprenons le fil interrompu par Dattier.
Nous avons une liste TV de 10 000 temps de vie produite par le code que j'ai donné plus haut (code qui ne contient ni exp, ni log).
On peut prendre prendre les quantiles qui divisent l'échantillon ordonné en 20 classes de 500.

Code:

def quantiles(liste,q) :
    l=len(liste)
    liste.sort()
    quant=[liste[((k+1)*l)//q-1] for k in range(q)]
    return quant

ydat=[(k+1)/20 for k in range(20)]
tdat=quantiles(TV,20)
print(tdat)

[51, 102, 158, 218, 283, 357, 430, 509, 601, 697, 803, 920, 1057, 1219, 1403, 1628, 1895, 2287, 2937, 9992]

La lectrice attentive aura remarqué qu'ici la médiane est à 697 au lieu du 689 annoncé plus haut. Explication : j'ai relancé ma simulation.

À partir de ces quantiles, on peut "fitter" avec une fonction de répartition exponentielle de type f(t) = a + b*exp(-c*t). Comme f doit être une fonction de répartition, on doit avoir f(0) = 0, donc b = -a et la limite en plus l'infini de f doit être 1, donc a = 1. Finalement f(t) = 1 - exp(-c*t) et on a un seul paramètre c à optimiser.
Allons-y, et affichons pour voir le "fit".

Code:

def func(t,c) :
    return 1-np.exp(-c*t)

popt, pcov = curve_fit(func, tdat, ydat, p0=[1/100])

t=np.linspace(0,10000,100)
plt.plot(tdat, ydat, 'b.', label='données')
plt.plot(t,func(t,popt[0]),'r-', label="fit : 1 - exp(-{:.5f}*t)".format(popt[0]))
plt.xlabel('t')
plt.legend()
plt.show()

On retrouve comme paramètre c du "fit" le 1/1000. Rien de surprenant, n'est-ce pas ?

Oui, c'est pas mal, mais un observateur inconnu dirait probablement "C'est du grand n'importe quoi !".
Où avez-vous trouvé cette fonction "fit" ? Moi, mes fonctions, e les écris moi-même.
Mais enfin c'est pas mal quand-même.

Je remarque tout de même au passage qu'il y a deux points à éclaircir ou justifier :
1- vous déclarez une variable 'p' que vous appelez "probabilité", c'est la probabilité de quoi ? Je vous rappelle que l'énoncé précisait que 'je' n'avait aucune idée de quoi que ce soit, 'je' cherchait juste un modèle pour modéliser une situation réelle : la façon de numériser le phénomène de durée de vie sans usure.
2- vous décidez de modéliser ce phénomène avec une fonction exponentielle, d'où sortez-vous cette information ? Avez-vous testé d'autres fonctions ? J'imagine que c'est le rapport médiane/moyenne qui est très voisin de ln(2), je suis sûr que "je" qui fait cette recherche ignore ce détail. Mais comme vous êtes un matheux expérimenté, ce rapport vous a mis la puce à l'oreille.

Je rappelle que le seul problème dans ce fil est de simuler la loi exponentielle sans exp ni log. Autrement dit, de produire des échantillons de temps de vie obéissant à une loi exponentielle, sans utiliser exp ou log dans le code.
Votre question 2 est donc hors-sujet. Je suppose que votre "je" n'est pas complètement idiot et a vu la décroissance exponentielle au lycée.
Pour la question 1

vous déclarez une variable 'p' que vous appelez "probabilité", c'est la probabilité de quoi ? Je vous rappelle que l'énoncé précisait que 'je' n'avait aucune idée de quoi que ce soit, 'je' cherchait juste un modèle pour modéliser une situation réelle : la façon de numériser le phénomène de durée de vie sans usure.

Je vous conseille de mieux lire ce que j'écris : "on choisit une probabilité p de disparaître par unité de temps (on prendra p petit)". Je répète : p est la probabilité de disparaître au cours de la prochaine unité de temps. Et comme vous voulez modéliser un phénomène sans usure, cette probabilité p est constante au cours du temps : pas d'usure qui augmenterait la probabilité de disparaître dans la prochaine unité de temps. Pourquoi p petit ? Parce qu'on veut modéliser un phénomène continu, on discrétise donc avec des petits pas.
Reste le choix de l'unité de temps. On suppose être informé de la demi-vie pour les durées de vie qu'on simule. On sait (parce qu'on a appris ses leçons au lycée) qu'avec notre modèle de simulation (probabilité p de disparaître par unité de temps), le nombre t d'unités de temps pour atteindre la demi-vie doit vérifier (1 - p)^t = 1/2, soit  t = - ln(2)/ln(1-p) ; comme p est pris petit, ln(1-p) est
très proche de -p et on peut prendre t = ln(2)/p.  Avec p = 1/1000 par exemple, ceci permet d'étalonner l'unité de temps en fonction de la demi-vie : l'unité de temps est la demi-vie divisée par 1000*ln(2).
Trève de baratin. Si vous avez un code pour produire des échantillons de durée de vie sans exp ni log, vous êtes bienvenu pour le communiquer.

Bonjour,

Code:

int main203() // suivant 203
{
// Simulation d'expérience exponentielle
  randomize();
  int ResExp[300];
  for (int i=0; i<300; i++)
  {
    int compt;
    for (compt=0 ; compt<2000; compt++)
    {
      int r=rand();
      if (r%79 == 0) break;
    }
    ResExp[i]=compt;
  }
  AfficheNormale(espion, ResExp, 300, "\nSimulation de loi exponentielle \n");
// on trie les valeurs, sens croissant
  TriSimple(ResExp, 300);
// on prend la moyenne par paquet de 15
  for (int i=0; i<300; i+=15)
  {
    float Moy=0.0;
    for (int j=i; j<15+i; j++)
    {
      Moy+=(float)ResExp[j];
    }
    Moy/=15.0;
    fprintf(espion,"%d  %0.3f\n",i+7,Moy);
  }
  fprintf(espion,"9999 9999\n0\n1\n9999\n");
  fclose(espion);
  return 0;
}



L'idée est que l'expérimentateur n'a aucune idée du résultat. Il ignore les notions de "probabilité" de "demi-vie". Donc il ne peut pas tricher comme vous le faites.
Comme il ne sait pas à quoi ça ressemble, il trie les résultats et simplifie en faisant des paquets de même taille auxquels il donne un numéro croissant suivant le rang.
Ensuite il fait examiner cette liste (rang moyenne) par son outil de calcul de régression, lequel choisi tout seul la "meilleure" formule. Cet outil sait le faire, parce que je lui ai appris.

Merci pour le code.
Les lignes

Code:

{
    int compt;
    for (compt=0 ; compt<2000; compt++)
    {
      int r=rand();
      if (r%79 == 0) break;
    }
    ResExp[i]=compt;
  }

montrent que votre échantillon de 300 durées de vie est donné par une loi géométrique de paramètre 1/79. Ce n'est donc pas différent de la manière dont je produis mon échantillon, voir ce qu'a rappelé Dattier ici : https://dlz9.forumactif.com/t2157-simuler-la-loi-exponentielle-sans-fonction-exp-ou-log#25534.. Seul varie le paramètre de la loi géométrique.
Ensuite, pour trouver le fit sur les données, vous travaillez avec les moyennes des tranches de 15. Pourquoi pas avec les quantiles ? Ça me semble pourtant plus logique, en termes de fonction de répartition : on compte le nombre de "morts" au temps t, c'est ce que font les quantiles.
Enfin, vous dites que votre outil choisit tout seul la "meilleure" formule. C'est très bien, mais c'est la meilleure parmi les quelques possibilités que vous avez fait figurer. Est-ce que votre outil essaie par exemple de fitter avec la formule
Comme vous pouvez le constater, ça marche pas mal :

Bonjour Dlzlogic.
Je rappelle ce que vous avez écrit :

Dlzlogic a écrit:Je rappelle que vous n'avez pas réussi à simuler une expérience satisfaisant cette loi [la loi exponentielle], sans utiliser vos connaissances mathématiciennes, c'est à dire simuler une loi géométrique qui approxime un loi exponentielle. Bref.

Que constate-t-on dans votre code ? Que vous simulez une expérience satisfaisant la loi exponentielle en simulant une loi géométrique de paramètre 1/79. Bref ...

Je vais vous répondre très simplement : 'je' et moi aussi ne savons pas ce qu'est une loi géométrique et ce qu'est un paramètre.
'Je' cherche seulement à simuler le phénomène. M. Jourdain ne savait pas qu'il écrivait en prose, mais il savait très bien ce qu'il voulait faire.

S'il vous plait, ne baratinez pas en niant l'évidence : vous répétez une épreuve de Bernoulli de probabilité de succès p = 1/79 et vous comptez le nombre d'épreuves pour aboutir à un succès. Ce nombre d'épreuves s'appelle une variable aléatoire de loi géométrique de paramètre 1/79.
Pourquoi me reprochez-vous alors de simuler une loi exponentielle en répétant une épreuve de Bernoulli de probabilité de succès p (p=1/1000 dans ce fil) et en comptant le nombre d'épreuves pour arriver au succès, alors que c'est exactement ce que vous faites ? Deux explications possibles :
1) vous ne savez pas lire,
2) vous êtes de mauvaise foi.
Je penche pour la première, puisque vous avez tout de même fini par donner votre code.

Ce nombre d'épreuves s'appelle une variable aléatoire de loi géométrique de paramètre 1/79.

Je ne connais ces termes et expressions que depuis que je lis des cours et des questions sur les forums. En particulier, ces distinctions entre discret et continue n'ont à mon avis aucune justification.
J'ai très bien compris que vous avez enseigné des choses suivant un certain programme et conformément à e que vous avez cru comprendre de certains documents, alors restez, s'il vous plait, à votre place.

Je ne connais ces termes et expressions que depuis que je lis des cours et des questions sur les forums.

Ce qui prouve bien que votre formation en probas a été très sommaire.
Alors, arrêtez de baratiner s'il vous plait. Des faits vérifiables, des codes, des nombres. Sur cette base, on peut argumenter.

Ben oui, vous continuez à mélanger "Théorie des partitions" et "Théorie des probabilités".

Encore et toujours du baratin. C'est pénible !
Oseriez vous prétendre que votre simulation ne consiste pas à répéter une épreuve de Bernoulli de proba de succès 1/79 et à compter le nombre d'épreuves pour arriver au premier succès ?

D'abord je ne sais pas ce qu'est une épreuve de Bernoulli. Si j'étais prof, je saurais l'expliquer à des élèves, mais je ne suis pas prof et ce qu'on appelle "épreuve de Bernoulli" est notion mathématique abstraite pour représenter un phénomène concret, ce qui est absurde.
Quand on veut simuler un phénomène avec un ordinateur, il n'y pas d'autre méthode que le comptage. Ce n'est pas pour autant que c'est une loi géométrique.

Une épreuve de Bernoulli est une épreuve avec deux issues possibles : succès avec probabilité p, échec avec probabilité 1-p.
Vous êtes en fait complètement ridicule en prétendant ne pas savoir ce qu'est une épreuve de Bernoulli  :
rand()%79 == 0 retourne vrai avec proba 1/79, , faux avec proba 78/79.

Vous n'avez toujours pas répondu à ma question

Pourquoi me reprochez-vous alors de simuler une loi exponentielle en répétant une épreuve de Bernoulli de probabilité de succès p (p=1/1000 dans ce fil) et en comptant le nombre d'épreuves pour arriver au succès, alors que c'est exactement ce que vous faites ?

Je vous rappelle mon code qui produit une liste de n temps de vie obtenus en comptant le nombre d'épreuves de probabilité de succès p pour arriver au premier succès :

Code:

def liste_temps_vie(p,n) :
    L=[]
    for _ in range(n) :
        t=0 ; vie=0
        while vie==0 :
            t+=1
            if rd.random() < p :
                vie=1
        L.append(t)
    return L

La commande

Code:

L=liste_temps_vie(1/79,300)
L.sort()
print(L)

produit :
[1, 1, 1, 2, 2, 2, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 14, 14, 14, 14, 15, 15, 16, 16, 16, 17, 18, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 21, 21, 21, 21, 23, 24, 24, 24, 25, 25, 25, 25, 27, 27, 27, 27, 27, 27, 27, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 37, 37, 38, 38, 39, 39, 39, 39, 39, 40, 40, 41, 41, 41, 41, 42, 42, 42, 44, 44, 44, 44, 45, 45, 45, 45, 46, 47, 47, 48, 48, 48, 48, 48, 50, 52, 52, 52, 52, 52, 55, 55, 56, 60, 61, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 64, 64, 65, 66, 67, 67, 67, 68, 69, 70, 71, 72, 72, 72, 72, 73, 74, 74, 75, 75, 76, 77, 78, 79, 79, 82, 83, 84, 85, 87, 87, 89, 89, 89, 90, 91, 92, 92, 93, 93, 93, 93, 96, 96, 98, 99, 100, 101, 101, 101, 102, 103, 103, 103, 103, 105, 107, 109, 109, 110, 112, 113, 114, 116, 116, 117, 117, 118, 119, 119, 121, 121, 124, 124, 125, 126, 127, 129, 130, 131, 132, 132, 133, 134, 135, 138, 145, 151, 152, 153, 154, 158, 159, 162, 162, 163, 167, 170, 173, 177, 185, 187, 190, 193, 194, 198, 206, 206, 210, 210, 211, 214, 218, 219, 219, 221, 224, 232, 235, 243, 247, 252, 282, 287, 289, 299, 305, 307, 315, 343, 351, 353, 356, 359, 407, 470, 501, 518, 636]
Les moyennes par tranches de 15 sont
[3.47, 9.73, 15.0, 19.87, 25.93, 30.8, 35.07, 38.93, 44.2, 51.73, 63.07, 70.87, 82.6, 95.2, 107.33, 122.2, 144.6, 185.4, 234.27, 387.27]
et les quantiles pour une division en 20 classes sont
[6, 12, 18, 24, 28, 34, 36, 41, 48, 61, 67, 75, 90, 101, 116, 131, 162, 210, 289].
Le QQ-plot par rapport à la loi exponentielle pour ces quantiles :


​

Il semble que devrais détailler un peu mon expression "je ne sais pas ce qu'est une variable de Bernoulli".
Pour cela, il me parait utile d'expliquer un peu ce que veut dire probabilités.
D'abord la définition : le rapport du nombre de cas favorables sur le nombre de cas possibles. Il est clair que rien ne dit que "nombre" doit être un nombre entier. Cela peut être un nombre réel, par exemple une aire, une mesure, c'est à dire quelque-chose qui peut être caractérisé par un nombre.
L'un de ces cas, pris indépendamment des autres n'a que très peu de signification, en d'autres termes il ne présent aucun intérêt. Il me semble inutile de la caractériser.
Il peut effectivement être intéressant d'étudier la probabilité du moment d'arrivée du "premier succès". Lorsque toutes les épreuves sont indépendantes, surtout si elles portent sur des objets différents, comme le contenu du porte-feuille du nième client dans une file d'attente ou la nième ampoule ou le nième atome de carbone 14, sauf pour des études particulières, c'est un aspect particulier de l'utilité des probabilités.
Par contre, les utilisations importantes des probabilités résident dans l'examen de choses provenant de la même source. Le contrôle de fabrication ou la statistique par sondage sont des exemples courants. Isoler les cas ne présente aucun intérêt, sauf éventuellement dans le cas d'exercices scolaires.

Quel baratineur vous faîtes !

Je ne dois pas être très bon, puisque vous ne comprenez rien de ce que j'explique.

La seule chose concrète que vous avez apporté dans ce fil est votre code de simulation. Rassurez-vous, je l'ai très bien compris.
Pour le reste, ça ne m'intéresse absolument pas de polémiquer sur votre baratin, qui tourne en boucle autour de "vous ne comprenez rien".