Shared5.sagewsOpen in CoCalc

Thème 5 : Simulation de variables aléatoires, estimation

Sujet de devoir : voir à la fin du document

Documentation Sagemath :

La documentation Sage en ligne

Quick Reference Statistics

Histogrames avec Sagemath

Premiers documents :

Histogramme pour la loi binomiale et graphe de la densité gaussienne avec SageMath (document pour le cours de proba en L2Mass)

TP Scilab de 2006 : TP0,2,3,4. Suivre le lien vers le cours pour le sujet des TP.

Le TP sur la méthode de Montecarlo en 2015-16 (scilab)

0. R (cf cours "Info - R" en LMASS) dans Sagemath

%r
#Avec R mais semble t-il pas d'échange entre R et Sage malgré ce qu'indique http://doc.sagemath.org/html/en/reference/interfaces/sage/interfaces/r.html
x=runif(10,0,1);x;summary(x)
  1. 0.187417119042948
  2. 0.464449916733429
  3. 0.974899226799607
  4. 0.928376206662506
  5. 0.179594443412498
  6. 0.60894227353856
  7. 0.65537861129269
  8. 0.116265019169077
  9. 0.197909407550469
  10. 0.711789044318721
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.1163 0.1900 0.5367 0.5025 0.6977 0.9749
%default_mode sage # change r to any mode
x

x
%r
x
  1. 0.187417119042948
  2. 0.464449916733429
  3. 0.974899226799607
  4. 0.928376206662506
  5. 0.179594443412498
  6. 0.60894227353856
  7. 0.65537861129269
  8. 0.116265019169077
  9. 0.197909407550469
  10. 0.711789044318721
L = r.rbinom(100,10,.5)
#ne marche pas
#cf https://ask.sagemath.org/question/8848/how-do-i-generate-a-random-number-according-to-the-binomial-distribution/

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 982, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> AttributeError: 'function' object has no attribute 'rbinom'

1. Stat descriptives avec Sagemath

#Stat descriptive, cf sage.stats.basic_stats? ou http://doc.sagemath.org/html/en/reference/stats/sage/stats/basic_stats.html
#sage.stats.basic_stats?
l=[1,1,2,3,3,4,5]
mode(l)
mean(l),mean(l).n(digits=3)
std(l), std(l).n(digits=3)
median(l)
P=frequency_distribution(l).function();P;P[1]
[1, 3] (19/7, 2.71) (sqrt(47/21), 1.50) 3 {1: 0.285714285714286, 2: 0.142857142857143, 3: 0.285714285714286, 4: 0.142857142857143, 5: 0.142857142857143} 0.285714285714286

2. Lois de proba discrètes ou continues. Cf :

http://doc.sagemath.org/html/en/prep/Quickstarts/Statistics-and-Distributions.html#distributions

http://doc.sagemath.org/html/en/reference/probability/sage/probability/probability_distribution.html

très technique http://doc.sagemath.org/html/en/reference/misc/sage/misc/randstate.html

http://doc.sagemath.org/html/en/reference/misc/sage/misc/prandom.html

#Simulation loi uniforme sur un intervalle
#avec la mention 'seed=1' les execution de cette cellules donneront toujours le même résultat, cf hps//ask.mathor/quetion/10194/seed---functions/
unif=RealDistribution('uniform', [0,2],seed=1)
[unif.get_random_element() for i in range(10)]
#Egalement uniform?
[0.834043996874243, 1.99436961626634, 1.440648978576064, 1.8651147224009037, 0.0002287621609866619, 0.25624889554455876, 0.60466513549909, 1.9980810307897627, 0.2935117851011455, 0.4721779525279999]
#cf random?
[2*random()+1 for i in range(10)]
[1.9087383406327822, 1.9811235568612375, 2.233497636857148, 1.3199847758534806, 1.6884498332296403, 1.691359264651321, 2.861548602342093, 1.7475436745236606, 2.7446343229519354, 1.6086455876789922]
#Simulation de la loi uniforme sur les n premiers entiers
#Cf ZZ.random_element?
#https://ask.sagemath.org/question/23377/fast-uniform-random-integer-generation/

Simulation d'une loi qcq sur l'ensemble \{0,\ldots,n\}

P = [0.1, 0.6, 0.3]
X = GeneralDiscreteDistribution(P)
l=[X.get_random_element() for i in range(100000)];l[:50]

[2, 0, 2, 1, 2, 1, 1, 2, 1, 2, 0, 1, 1, 0, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 0, 2, 1, 0, 1, 0, 2, 2, 2, 0, 2, 1, 2]
#Fréquence de 0,1,2 dans l
for k in [50,1000,100000]:
    print k,": ",[sum(l[i]==j for i in range(k))/k.n(digits=3) for j in [0,1,2]]
50 : [0.140, 0.480, 0.380] 1000 : [0.0930, 0.618, 0.289] 100000 : [0.0987, 0.602, 0.299]
#Représentation graphique de la loi discrète P, des fréquences de l : diagramme en batons
from sage.plot.histogram import Histogram
bar_chart([sum(l[i]==j for i in range(20))/20 for j in [0,1,2]],width=0.2,rgbcolor=(0,1,0))+bar_chart(P, width=0.1)

Simulation d'une va de Bernouilli, intervalle de confiance de la fréquence observée de 1

Que rend chacune des instructions suivantes ? Comment calculer la proportion de 1 dans la liste rendue par la dernière instruction ?

(random()<.3).real
p=.3;[(random()<p).real for _ in range(100)]
0 [0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
#La 1ère instruction donne 1 si random() a donné un nbre inférieur à 0.3, 0 sinon (donc 1 avec proba 0.3), la 2ème donne une liste de 100 résultats d'exécution de la 1ère instruction
l=[(random()<p).real for _ in range(100)];f=Integer(sum(l))/Integer(len(l));l;f;f.n(digits=2) #sum(l)/len(l) est de type int ; la division de deux objets de type int est de type int (Python 2.7), voir ci-dessous 
[0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0] 17/50 0.34
int(1)/int(2)
len([1])/len([1,2])
type(len([1]))
Integer(len([1]))/Integer(len([1,2]))
type(Integer(len([1])))
0 0 <type 'int'> 1/2 <type 'sage.rings.integer.Integer'>

Soit p\in[0,1] et (x_k)_{1\leq k\leq N} une suite d'éléments de \{0,1\} obtenue par simulation d'une va de loi \mathcal{B}(p) . Notons f_1 la proportion (ou fréquence) de 1 observée dans la suite (x_k) ; N*f_1 est une variable aléatoire de loi binomiale \mathcal{B}(N,p) ; f_1 suit une loi proche de la loi normale \mathcal{N}(p,\sqrt{p(1-p)}/\sqrt{N}) , Cf le théorème de Moivre-Laplace, cas particulier du Théorème central limite.

On connaît les quantiles de la loi normale : si Y est de loi \mathcal{N}(\mu,\sigma) alors par exemple \mathrm{P}(X\in[\mu-1.96\sigma, \mu+1.96\sigma])\simeq 0.95

En utilisant la feuille 1, pouvez vous vérifier cette affirmation ?

On s'attend donc à ce que, si on répète le calcul de f_1 un certain nombre de fois, 95% des valeurs pour f_1 obtenues soient dans l'intervalle

[p-1.96\frac{\sqrt{p(1-p)}}{\sqrt{N}},\ p+1.96\frac{\sqrt{p(1-p)}}{\sqrt{N}}]

Prenons p=.15, N=100 ; pouvez vous vérifier expérimentalement cette affirmation ? Voir le sujet et corrigé du TP3 de 2006

var('m s')
f(x)=exp(-x^2/2)/sqrt(2*pi);f #densité de la loi normale centrée réduite
g(m,s,x)=f((x-m)/s)/s #densité de la loi normale d'espérance m d'écart type s
integrate(g(m,s,x),x,m-1.96*s,m+1.96*s).n()
print
p=.15;n=100;l=[frequency_distribution([(random()<p).real for _ in range(n)]).function()[1] for _ in range(200)];[i.n(digits=2) for i in l]
s=sqrt(p*(1-p)/n);s
frequency_distribution([i>p-1.96*s and p<p+1.96*s for i in l]).function()[true]
(m, s) x |--> 1/2*sqrt(2)*e^(-1/2*x^2)/sqrt(pi) 0.950004209703559 [0.23, 0.15, 0.22, 0.19, 0.15, 0.16, 0.18, 0.16, 0.17, 0.19, 0.060, 0.10, 0.12, 0.21, 0.22, 0.19, 0.13, 0.14, 0.12, 0.17, 0.23, 0.070, 0.12, 0.17, 0.13, 0.19, 0.13, 0.19, 0.24, 0.16, 0.15, 0.19, 0.14, 0.12, 0.16, 0.16, 0.16, 0.14, 0.17, 0.17, 0.090, 0.090, 0.18, 0.17, 0.14, 0.11, 0.090, 0.13, 0.18, 0.24, 0.19, 0.14, 0.13, 0.17, 0.16, 0.090, 0.19, 0.12, 0.16, 0.20, 0.15, 0.080, 0.090, 0.21, 0.13, 0.17, 0.17, 0.14, 0.17, 0.060, 0.20, 0.17, 0.14, 0.14, 0.18, 0.17, 0.16, 0.19, 0.20, 0.15, 0.15, 0.18, 0.17, 0.18, 0.13, 0.15, 0.18, 0.18, 0.15, 0.12, 0.12, 0.19, 0.080, 0.18, 0.15, 0.17, 0.13, 0.17, 0.17, 0.19, 0.15, 0.16, 0.16, 0.18, 0.15, 0.23, 0.20, 0.15, 0.090, 0.12, 0.14, 0.16, 0.080, 0.090, 0.11, 0.18, 0.16, 0.17, 0.16, 0.15, 0.19, 0.16, 0.15, 0.12, 0.19, 0.15, 0.11, 0.14, 0.23, 0.12, 0.13, 0.18, 0.19, 0.21, 0.19, 0.19, 0.17, 0.25, 0.14, 0.15, 0.12, 0.11, 0.11, 0.15, 0.23, 0.20, 0.14, 0.17, 0.12, 0.16, 0.20, 0.14, 0.080, 0.21, 0.18, 0.11, 0.13, 0.12, 0.13, 0.15, 0.20, 0.17, 0.19, 0.19, 0.19, 0.14, 0.13, 0.18, 0.22, 0.18, 0.12, 0.17, 0.19, 0.15, 0.17, 0.21, 0.20, 0.16, 0.15, 0.14, 0.17, 0.15, 0.21, 0.13, 0.12, 0.13, 0.18, 0.12, 0.22, 0.14, 0.11, 0.20, 0.14, 0.17, 0.16, 0.18, 0.12, 0.13, 0.14, 0.11] 0.0357071421427143 0.965000000000001

Devoir : Simulation d'une chaîne de Markov

Sur les chaînes de Markov, voir le cours de l'option "Markov", ou cette page (cours en 2013-2016).

Ecrire une fonction qui prend comme argument une matrice stochastique A de taille n (avec la convention que chaque colonne somme à 1 plutôt que chaque ligne), un élément x_0 de \{1,\ldots,n\} , un entier N et rend une liste (x_k)_{0\leq k\leq N} d'éléments de \{1,\ldots,n\} x_k est déduit de x_{k-1} et d'une expérience aléatoire telle que \mathrm{P}(x_k=i|x_{k-1}=j)=A_{i,j}

Ecrire une fonction qui prend comme argument une liste (x_k) d'éléments de \{1,\ldots,n\} et qui rend une estimation statistique de la matrice A dont (x_k) pourrait être déduit par la première fonction.

Représenter graphiquement la liste produite

Tester ces fonctions avec les matrices

A=\left(\begin{array}{ll} 0.5 & 0.5\\ 0.5& 0.5\end{array}\right),\qquad A=\left(\begin{array}{ll} 0.2 & 0.8\\ 0.8& 0.2\end{array}\right)

Qu'est ce qui distingue les représentations graphiques ?

Sujet au tableau