Réponse au TP de Statistique (projet)

1.1 rand() a produit successivement
  0.2113249, 0.7560439, 0.0002211, 0.3303271, 0.6653811, 0.6283918
rand() est une fonction ne prenant pas d'argument et rendant un résultat non déterministe.

1.2 a=1:1000;for i=1:1000,a(i)=rand(),end
On voit l'affichage de la valeur de a(1) puis de l'indice i de 2 à 50.
Premier graphique pour 30 et 200 itérations :


La suite de points remplit le segment [(0,0),(1,0)] de façon non uniforme en apparence mais tendant vers l'uniforme lorsque le nombre de points affichés grandit. La qualité de l'affichage ne permet pas de distinguer les graphiques correspondant à un grand nombre d'itérations.

Histogramme :



1.3
for i=1:1000,a(i)=2*rand()-1,end
On voit s'afficher :
-->for i=1:1000,a(i)=2*rand()-1,end
 a  =
 
 
         column 1 to 5
 
!   0.1215891    0.0683740    0.5608486    0.6623569    0.7263507 !
 
         column  6 to 10
 
!   0.1985144    0.5442573    0.2320748    0.2312237    0.2164633 !

Il y a manifestement un bug dans l'affichage du résultat, mais si on entre ensuite l'instruction a (et on valide), on obtient les 25 premiers termes de la suite a avec le résultat attendu : nombres compris entre -1 et 1 :
-->a
 a  =
 
 
         column 1 to 5
 
!   0.1215891  - 0.2742239  - 0.3823348    0.9626830    0.4856965 !
 
         column  6 to 10
 
! - 0.6767217  - 0.2694956  - 0.2329859  - 0.7661638    0.7822492 !

  b=1:1000;for i=1:1000,b(i)=2*rand()-1,end;
  xbasc();for i=1:1000,plot2d(a(i),b(i),0,rect=[-1.2,-1.2,1.2,1.2]),end
C'est un peu long. En fait la commande suivante fait la même chose mais beaucoup plus rapidement (chercher plot2d dans l'aide de Scilab).

  xbasc();plot2d(a,b,0,rect=[-1.2,-1.2,1.2,1.2])

Le résultat est satisfaisant. Si on redéfinit a par l'instruction
  for i=1:1000,a(i)=-1+2*rand()^2,end
on obtient

Les points s'accumulent d'avantage à gauche, la distribution n'est pas uniforme dans le carré.

2.1
bool2s(1<=2) donne 1 et bool2s(2<1) donne 0. On propose donc
  bool2s(x^2+y^2<=1)
qui rend une erreur parce que x (et y) non pas encore été définies.
x=sqrt(3)/2;y=1/2;bool2s(x^2+y^2<=1) donne 1 et x=sqrt(3)/2+0.001;y=1/2;bool2s(x^2+y^2<=1) donne 0 (on joue ici avec la précision de scilab dans le calcul avec les réels).

On redéfinit a comme en 1.3 puis c par
  c=1:1000;for i=1:1000,c(i)=bool2s(a(i)^2+b(i)^2<=1),end
sum(c) donne 792

Une approximation de pi est donnée par
   4*sum(c)/1000
On obtient 3.168

2.3 On recommence de façon synthétique par les deux instructions
  for i=1:1000,c(i)=bool2s(rand()^2+rand()^2<=1),end
  4*sum(c)/1000
 L'affichage est intempestif. On est plus heureux avec l'instruction
  c=bool2s(rand(1:1000)^2+rand(1:1000)^2<=1);4*sum(c)/1000
On obtient successivement
3.232, 3.076, 3.048, 3.068
 Les approximations ne sont pas très bonnes.

3.
d=1:400;for i=1:400,for j=1:1000,c(j)=bool2s(rand()^2+rand()^2<=1),end,d(i)=4*sum(c)/1000,end
mean(d) donne 3.14488. Pour juger de cette approximation, il faut executer plusieur fois l'instruction. On obtient successivement
3.14146, 3.14782, 3.13919, 3.13713
Ce qui semble donner une approximation de pi à 0.01 près.

4.1
x=2.97:0.02:3.31;de=d(2.97<=d & d<=3.31);xbasc();histplot(x,de)

4.2
m=%pi
s=sqrt(m*(4-m)/1000)
function y=f(x);y=1/(s*sqrt(2*%pi))*exp(-.5*((x-m)/s)^2);
endfunction
fplot2d(2.97:0.001:3.31,f)

Le graphe de f épouse assez bien le bord de l'histogramme.

5.1 se=sqrt(variance(d)) donne 0.0497505 qu'on compare avec la valeur idéale s = 0.0519304

5.2
e=1:400;for i=1:400,e(i)=bool2s(d(i)-se<=%pi & %pi<=d(i)+se),end
sum(e) donne 271

Le pourcentage cherché est donné par sum(e)/400 ce qui donne 0.6775 soit environ 68%.

5.3 avec 2se on obtient 0.95 soit 95%
5.4
function y=g()
for i=1:400,c=bool2s(rand(1:1000)^2+rand(1:1000)^2<=1),d(i)=4*sum(c)/1000,end
se=sqrt(variance(d))
for i=1:400,e(i)=bool2s(d(i)-se<=%pi & %pi<=d(i)+se),end
y=sum(e)/400
endfunction

 Les exécussions succéssives de l'instruction g() donnent
0.6825, 0.7075, 0.695, 0.6575
etc.

F-X Dehon, 21 décembre 2005