previous up next_inactive
Up: Return to previous menu

Ensait - E2 - Aide à la décision


Date: Corrigé de l'évaluation du 2008-01-18

durée 2h, tous documents autorisés

A nouveau, quelques remarques d'ensemble

  1. La première qualité que l'on attend d'un ingénieur est de savoir présenter ses conclusions. Les listings, les calculs, les graphiques et autres "sorties d'ordinateur" ne peuvent en aucun cas remplacer un relevé de conclusions, rédigé de façon précise et scientifique.
  2. Il est donc nécessaire d'ordonner différents types de matériaux de façon à produire un document organisé, produisant une liste de réponses argumentées à la liste des questions de l'énoncé.
  3. La méthode RECOMMANDÉE était de rédiger les réponses à la main et de faire des couper-coller avec des ciseaux et de la colle pour insérer les "sorties d'ordinateur". L'expérience confirme que cette méthode est à la fois la plus rapide et la plus robuste. C'est aussi celle qui permet de se concentrer sur les compétences évaluées et non sur les compétences/incompétences en dactylographie.
  4. Principe d'adéquation : lorsque la question change, il est assez fréquent que la réponse adéquate change aussi. Quand on constate que le graphe théorique du sujet de révision n'est pas en accord avec le graphe expérimental de l'évaluation, ce n'est pas la faute à Walker. C'est que les données ont changé.
  5. Un minimum de rigueur orthographique et grammaticale est attendu. La "varaince", "on génère les nombre 1, 2, autre", "supperposer", "le nombre d'unité de temps avant l'arrivé du", "algotrythme" et le "test du Qui2" créent une certaine surprise à un niveau master. Quant à "on résonne de la même manière"...
  6. Enfin, les copies sont évaluées en proportion de leur contenu et pas en fonction de leur poids. Dans cet ordre d'idées, réimprimer vingt pages de bibliothèques (dont la simulation de l'ensemble des files d'attente d'une agence bancaire) n'est pas utile.
  7. Les notes s'étagent de 03 à 20, avec une moyenne de 11.5. Huit étudiants ont 17 ou plus. Sept étudiants ont 07 ou moins.

Rappel des consignes données dans le sujet

  1. Les valeurs initiales ont changé par rapport à des exercices précédents. Savoir où et comment modifier les programmes utilisés fait partie des compétences évaluées.
  2. Les remarques faites à l'occasion du module "Plans d'expérience" continuent d'être valables. Elles sont consultables à l'adresse  http://www.douillet.info/~douillet/cours/planx_ds07a/index.html.
  3. Par ailleurs, l'attention des étudiants est attirée sur le fait que le trafic réseau de leur ordinateur est susceptible d'être enregistré pendant la durée de l'évaluation.
  4. Les bibliothèques decis et waits doivent être à jour, c'est à dire avoir été téléchargées récemment sur le site web. Les numéros de version sont :

    $\displaystyle decis\ge04 ; waits\ge18$

  5. Modifier la ligne $ 8$ de chacun de ces deux programmes pour qu'elle indique votre répertoire de travail.
  6. Avant chaque simulation, réinitialiser le générateur aléatoire par la commande :
    grand('setsd',2ddmmyy)
    où votre date de naissance est dd/mm/19yy

1 Algorithme de Walker

On s'intéresse au générateur à quatre choix décrit dans decis.sce. Initialiser le générateur aléatoire. Puis utiliser initpick pour introduire les valeurs $ N=500$, $ p_{1}=0.15$, $ p_{2}=0.4$ et $ p_{3}=0.25$ (vérifier les valeurs introduites en examinant le message affiché dans scilex).
  1. Exécuter pickc (le générateur de Walker). Que valent alors les vecteurs gap et aux ?
    1. Ne pas oublier d'initialiser le générateur 
      grand('setsd',123456);
    2. initpick(); pr=[p1,p2,p3,p4]; pickc(); 
      printf( 'gap= '); printf('%5.3f ',gap'); 
      printf('\nalt= '); printf('%5d ',alt'); printf('\n\n');
    3. On obtient : 
      N=500, p1=0.15, p2=0.40, p3=0.25, p4=0.20, seed=123456  
      gap= 0.600 1.000 1.000 0.800  
      alt= 2     2     3     2     .
  2. Expliquer comment fonctionne le générateur gen, et montrer qu'il engendre effectivement les nombres $ 1,2,3,4$ avec les $ p_{1}, p_{2}, p_{3}, p_{4}$ comme fréquences théoriques. On pourra utiliser sci2exp(gen) pour disposer d'un affichage commode de la fonction gen.
    1. printf('%s\n', sci2exp(gen)) 
      qui donne : 
      createfun(["y=%fun";  
      " z = rand(); y = ceil(z * 4)";  
      " if (y - 4 * z) > gap(y) then y = alt(y); end";  
      ])
    2. Au niveau bac+4, expliquer comment fonctionne le générateur gen ne consiste pas à traduire le code en "prendre rand, puis multiplier par 4, puis prendre la partie entière, puis...". Cela consiste à expliquer pourquoi faire tout cela va effectivement conduire aux fréquences voulues. La base de ce fonctionnement est que :
      z est une variable aléatoire réelle, uniformément distribuée sur $ \left[0,\,1\right]$
      c$ \doteq$ceil(z*4) est une variable entière uniformément distribuée sur $ \left\{ 1,\,2,\,3,\,4\right\} $
      v$ \doteq$ceil(z-4)-4*z est une variable aléatoire réelle, uniformément distribuée sur $ \left[0,\,1\right]$
    3. Dans le cours, c est le candidat, v est le verdict.
    4. On obtient $ Pr\left( Y=y \right) $ par la formule des probabilités totales :

      $\displaystyle Pr\left( Y=y \right) =\sum_{c=1}^{c=4}Pr\left( Y=y\mid C=c \right) \times Pr\left( C=c \right) $

      Les $ Pr\left( C=c \right) $ valent toutes $ 1/4$ tandis que, par définition,

      $\displaystyle Pr\left( Y=c\mid C=c \right) =gap(c),\quad Pr\left( Y=alt\left(c\right)\mid C=c \right) =1-gap(c)$

    5. chances=gap; 
      for j=1:4 do chances(alt(j))=chances(alt(j))+1-gap(j); end 
      printf('\nchances= '); printf('%5.3f ',chances'); 
      printf('\n pr= '); printf('%5.3f ',pr'); printf('\n'); 
      donne :
      chances= 0.600 1.600 1.000 0.800  
      1]6.5mm 
        pr= 0.150 0.400 0.250 0.200 
      prouvant que les chances de sélection d'un candidat sont bien celles qui conviennent.
  3. Le programme pickc trace automatiquement l'histogramme des $ N=500$ valeurs engendrées. Ajouter le tracé de l'histogramme des fréquences théoriques.
    1. Utiliser histplot(bor, [1,2,2,3,3,3,4,4,4,4]) pour tracer l'histogramme théorique était légèrement inapproprié : cette "ruse" était adaptée à l'exercice de préparation, où les probabilités étaient : $ 0.1,\,0.2,\,0.3,\,0.4$.
    2. On pouvait, par contre, utiliser une liste à 20 valeurs : 
      histplot(bor, [1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,4,4,4,4])
    3. Et même (génération, puis test) :
      li=[1*ones(1:3),2*ones(1:8),3*ones(1:5),4*ones(1:4)]  
      tmp=tabul(li,'i'); tmp=tmp(:,2)'; tmp/sum(tmp)  
      histplot(bor,li)
    4. Comme les deux histogrammes sont faits pour se ressembler, il n'était pas inutile de passer l'histogramme expérimental en "gras".
      curax=gca(); curax.children(2).children.thickness=2;
    Figure 1: Algorithme de Walker
    \includegraphics[height=60mm,keepaspectratio]{figures/walker}

  4. La variable lesx contient les valeurs engendrées. Déterminer combien de fois les valeurs $ 1, 2, 3, 4$ ont été obtenues expérimentalement. Utiliser un test du $ \chi^{2}$ pour comparer avec les valeurs théoriques.
    1. On pouvait utiliser
      tabul(lesx, 'i'), obs=ans(:,2),hyp=pr*N
    2. On pouvait utiliser aussi
      [ind,obs]=dsearch(lesx,bor); hyp=pr*N; 
      printf('\nvaleur'); printf('%5d ',(1:4)')  
      printf('\n hyp'); printf('%5d ',round(hyp'))  
      printf('\n obs'); printf('%5d ',obs'); printf('\n');
    3. Par les deux méthodes, on obtient :
      valeur  1  2  3  4  
      1]7mm 
      hyp 75 200 125 100  
      1]7mm 
      obs 63 212 117 108
    4. L'indicateur de Pearson est donné par
      pearson=sum((hyp-obs).^2 ./ hyp); 
      Remarque : pour comprendre l'effet de l'opérateur "./", il peut être utile d'écrire la commande précédente sous la forme :
      (hyp-obs).^2 ./ hyp, pearson=sum(ans) 
      En tout cas diviser par hyp (en utilisant l'opérateur "/") ne conduit pas au résultat: le logiciel ne veut pas croire que vous avez été suffisamment inattentif pour vouloir utiliser l'inverse d'une matrice de taille $ 1\times4$ et traduit cela en "pseudo-inverse"... concept n'ayant rien avoir avec le calcul à entreprendre.
    5. On peut alors calculer la variable réduite associée et se référer à la courbe donnée en cours. On peut aussi calculer la probabilité cumulée associée à la valeur de l'indicateur de Pearson. Dans tout les cas il est "inattendu" de traiter la valeur de $ \chi_{red}^{2}$ comme une variable de Gauss : la loi du chi-square n'est pas la loi normale.
    6. Les commandes :
      pearson_red=(pearson-3)/sqrt(6); prob=cdfchi('PQ',pearson,3);  
      printf('\n pearson= %5.3f, pearson_red= %5.3f, prob= %f\n',  
      1]5mm 
      pearson, pearson_red, prob);  
      conduisent à :
      pearson= 3.792, pearson_red= 0.323, prob= 0.715182
    7. Bilan : la fonction gen est un générateur acceptable (remarque : tel était le but poursuivi...).

2 Test de Student-Fisher

Remarque préalable : le métier d'ingénieur ne consiste pas à cliquer toujours et systématiquement sur le même mulot, et le test de Student-Fisher ne consiste pas à cliquer hors de propos sur le corrigé de l'exemple 3.3.6 http://www.douillet.info/~douillet/cours/decis/node4.html. On considère à nouveau la variable aléatoire $ X\in\left\{ 1, 2, 3, 4\right\} $ avec $ Pr\left( X=i \right) =p_{i}$.
  1. Déterminer l'espérance et la variance de $ X$.
    1. A nouveau, veut dire que l'on continue avec la même variable. Les probabilités continuent d'être les mêmes $ p_{i}$ que précédemment.
    2. On a le choix entre
      mu=mean(li); sigma=sqrt(mean(li.^2)- moy^2);
      et
      qui=[1,2,3,4]; 
      mu=meanf(qui,pr); sigma=sqrt(meanf(qui .^2,pr)-moy^2);
    3. Dans les deux cas, la commande variance ne donne pas la variance (cet état de fait est signalé dans le cours).
    4. printf('\nmu= %f, sigma^ 2=%f\n', mu, sigma^ 2);
      donne
      mu= 2.500000, sigma^ 2=0.950000
  2. Déterminer l'intervalle de confiance (aux risques de $ 5\%$ puis de $ 1\%$) de la moyenne correspondant à un échantillon de taille $ N=11$ puis $ N=500$.
    1. La formule voulue est :

      $\displaystyle moy\in\left[\mu\pm k\left(\alpha\right)\,\sqrt{\frac{\sigma^{2}}{N}}\right]\qquad;\qquad k\, Gauss$ (1)

      lorsque l'on encadre la moyenne de l'échantillon à partir des paramètres de la population et

      $\displaystyle \mu\in\left[moy\pm t\left(\alpha\right)\,\sqrt{\frac{var_{est}}{N}}\right]\qquad;\qquad t\, Student$ (2)

      lorsque l'on encadre la moyenne de la population à partir des paramètres de l'échantillon.
    2. Pour un encadrement symétrique, un risque de $ 5\%$ se répartit en $ 2.5\%$ de chaque côté. La probabilité à utiliser dans les cdf est donc $ p=0.975$. Pour le risque $ 1\%$, on utilise la probabilité cumulée $ q=0.995$.
    3. Pour ne pas tout écrire en double, on utilise une boucle 
      for n=[500,11] do for p=[0.975, 0.995] do 
      k=cdfnor('X',0,1,p,1-p); delta=k*sigma/sqrt(N);  
      printf('au risque %d%%, k= %4.2f et moy(%3d) in [%6.4f,%6.4f]\n', .. 
      1]5mm 
      round(200*(1-p)), k, N, mu-delta, mu+delta)  
      end; end
    4. Les valeurs numériques correspondant à (1) sont :
      au risque 5%, k= 1.96 et moy(500) in [2.4146,2.5854]  
      au risque 1%, k= 2.58 et moy(500) in [2.3877,2.6123] 
      au risque 5%, k= 1.96 et moy( 11) in [1.9240,3.0760] 
      au risque 1%, k= 2.58 et moy( 11) in [1.7430,3.2570]
  3. Comparer avec la moyenne de l'échantillon déjà obtenu ($ N=500$). Simuler en outre un échantillon de taille $ N=11$ et calculer sa moyenne.
    1. Pour l'encadrement expérience vers théorie, il faut utiliser
      t=cdft('T',n-1,p,1-p); delta=t*sqrt(var/n);
      et on obtient (pour $ N=500$)
      N=500 et risque 5%, t= 1.96 et mu in [2.4551,2.6249]  
      N=500 et risque 1%, t= 2.59 et mu in [2.4282,2.6518]
    2. Pour $ N=11$, nous sélectionnons lesx(10:20). Les calculs donnent : 
      N= 11 et risque 5%, t= 2.23 et mu in [1.6202,2.5616]  
      N= 11 et risque 1%, t= 3.17 et mu in [1.4214,2.7604]
    3. Le facteur de Student est (quasiment) égal au facteur normal lorsque $ N$ est grand. Par contre, pour $ N$ petit, $ t$ est plus grand que $ k$ de façon à tenir compte des risques de sous-estimation de l'écart-type de la population.

3 File d'attente

On utilise désormais le programme de simulation de files d'attentes décrit dans waits.sce.

3.1 Petit essai

  1. Modifier la fonction petitessai de la façon suivante: (ligne 6 à 8)
    grand('setsd',2ddmmyy)  
    deff('s=arv',sprintf('s=grand(1, ''uin'', %f, %f)',10, 19)) 
    deff('s=srv',sprintf('s=grand(1, ''uin'', %f, %f)',16, 26))
  2. Lancer le programmer et en commenter le résultat.
    Figure 2: Un petit essai
    \verbatiminput{figures/petitessai.txt}

    1. On obtient le Figure 2.
    2. On remarque que $ \rho=\left(16+26\right)/\left(10+19\right)\approx1.45$. Le système est donc congestif. Abandonné à lui même, la file d'attente va tendre à croître indéfiniment.
    3. Le principal commentaire à faire est de dater chaque événement. Le délai temporel entre deux événements correspond, à partir d'un événement donné, au plus proche des deux événements à venir qui est soit une arrivée, soit un départ :

      $\displaystyle delai=\min\left(nextarv,\, agenda(1)\right)$

    4. numcli est le numéro du client dernier arrivé. On obtient le numéro d'arrivée des autres clients en partant de la fin. charge est le temps de séjour résiduel du dernier client arrivé. Au point de vue du serveur, c'est la charge de travail résiduelle autrement dit le temps à passer avant d'avoir fini de servir tous les clients actuellement présents. nbsej n'est pas le nombre de clients en attente, mais le nombre total de clients présents.
    5. A la date $ 010$, le client $ \char93 2$ arrive pour un séjour $ \Delta t=37$. On en déduit que sa date de départ sera $ 10+37=47$ (ce qui est bien le cas). Ce séjour est la somme de l'attente et du service. Comme l'attente est liée au séjour résiduel du client précédent, et que celle-ci vaut $ 16$, on en déduit que le service est $ 21$. Ceci se retrouve comme étant le délai séparant le départ des clients $ \char93 1$ et $ \char93 2$ soit $ 47-26=21$.
    6. On peut même remarquer le lien avec "Fig 5.4: Graphical proof of the Little's theorem, from Sinclair, 2006", http://www.douillet.info/~douillet/cours/decis/node6.html

3.2 Simulation M/Ga/1

Utiliser le programme 'système M/Ga/1' avec les réglages $ marv=65, msrv=59, a=2.4, N=50000$ et $ seed=2ddmmyy$. Le processus d'interarrivée $ X$ est alors une loi exponentielle de paramètre $ \lambda$ et le processus de service $ Y$ est une loi Gamma de paramètres $ a$ et $ b$ tels que $ \mathrm{E}\left(X\right)=marv$ et $ \mathrm{E}\left(Y\right)=msrv$.
Figure 3: Les graphes engendrés par la simulation
\includegraphics[width=80mm]{figures/durees} \includegraphics[width=80mm]{figures/sejour}
\includegraphics[width=80mm]{figures/service} \includegraphics[width=80mm]{figures/resid}

  1. Déterminer le débit d'arrivée $ \lambda$ à partir de $ marv$ et donner la distribution (pdf) de probabilité associée. Quelle est la variance du processus d'inter-arrivée ?
    1. Par définition

      $\displaystyle \lambda=\frac{1}{marv}=\frac{1}{65}\approx0.01538$

    2. Pour une loi exponentielle, le pdf (probability density function) est :

      $\displaystyle Pr\left( A\in\left[a,\, a+\mathrm{d}a\right] \right) =\lambda\,\exp\left(-\lambda\, a\right)\,\mathrm{d}a$

    3. La variance est $ \left(1/\lambda\right)^{2}=4225$.
  2. Rappeler quelle est la distribution de probabilité associée à la loi Gamma de paramètres $ a$ et $ b$. Déterminer la valeur de $ b$ en fonction de $ a$ et de $ msrv$. Quelle est la variance des services ?
    1. Une loi gamma réduite de paramètre $ a$ est

      $\displaystyle Pr\left( S\in\left[s,\, s+\mathrm{d}s\right] \right) =\frac{1}{\Gamma\left(a\right)}\, s^{a-1}\,\exp\left(-s\right)\,\mathrm{d}s$

      tandis qu'une loi gamma de paramètres $ \left(a,\, b\right)$ est

      $\displaystyle Pr\left( S\in\left[s,\, s+\mathrm{d}s\right] \right) =\frac{1}{\Gamma\left(a\right)}\,\frac{s^{a-1}}{b^{a}}\,\exp\left(-s/b\right)\,\mathrm{d}s$

    2. On a $ msrv=a\, b$ et donc $ b=msrv/a=59/2.4\approx24.58$
    3. La variance d'une loi Gamma est $ a\, b^{2}=msrv^{2}/a\approx1450.42$
    4. aaa=msrv^2/vsrv, bbb=msrv/aaa
  3. Rappeler ce qu'est le paramètre $ \rho$ et donner sa valeur.
    1. Ce paramètre est le rapport des débits

      $\displaystyle \rho=\frac{\lambda}{\mu}=\frac{msrv}{marv}=\frac{59}{65}\approx0.9077$

  4. Rappeler ce que signifient les graphes "1-Durée des états" et "2-États lors de l'arrivée". Pourquoi sont-ils si ressemblants dans l'exemple étudié ?
    1. Le graphe 1 donne le temps passé par le système dans les différents états d'occupation (c'est le point de vue du serveur).
    2. Le graphe 2 donne la statistique du nombre de clients en cours de séjour juste avant l'arrivée d'un nouveau client.
    3. Le fait que ces deux graphes se ressemblent n'est pas une règle générale, mais un cas spécial causé par le fait que les arrivées sont exponentielles.
  5. Le programme affiche diverses moyennes, écart-types et auto-corrélations. Pour ce qui est des services, donner l'intervalle de confiance à $ 95\%$ concernant la moyenne. Même question pour la variance (on signale que le $ M^{4}$ du cours vaut $ M^{4}=3a\left(a+2\right)b^{4}$ pour une loi Gamma).
    1. Les résultats de la simulation sont
      elapsed time = 19.579000 - N= 50000  
      durées totale=1612343.704758, durées à vide=139841.687500  
      nom         moyenne écart_type auto_corr. 
      durees   7.29979658 7.08189504 
      autres   7.28934213 7.06893418 0.98602434 
      sejours  470.7728069 422.9218919 0.98689272 
      services  58.9049995 38.3610466 0.00024394 
      residuels 42.0926403 36.0259904 0.30053120
    2. Ils peuvent être récupérés par mv=drawq(rep='1s-gam3');
    3. tmp=msscanf(mv(5,:),'%s %f %f%f') 
      z_moy_serv=(tmp(2)-msrv)/sqrt(vsrv/numcli)
      donne $ -0.3944$, valeur largement acceptable pour une variable normale.
    4. Le cours donne $ var\left(s^{2}\right)=\frac{1}{n}\left(M^{4}-\frac{n-3}{n-1}\,\sigma^{4}\right)$. On a donc :
      var_s2=3*aaa*(aaa+2)*bbb^4 -vsrv^2  
      z_var_serv=(tmp(3)^2-vsrv)/sqrt(var_s2/numcli)
      ce qui donne $ 1.0872$, valeur également acceptable.
  6. Pour ce qui est des services résiduels, considérer le fait que la loi en masse associée à la loi $ Gamma\left(a,b\right)$ est à son tour une loi Gamma, de paramètres $ a+1, b$ et obtenir les intervalles de confiance correspondants.
    1. Pour la répartition en masse des services, on obtient les valeurs théoriques : 
      mmas=83.583333, vmas=2054.756944
    2. Pour le service résiduel on a mres=mmas/2 (formule du cours) et vres=vmas/3 + mmas^ 2/12 (formule non demandée) :
      mres=41.791667, vres=1267.100116, sres=35.596350
    3. On compare avec les valeurs expérimentales : 
      residuels 42.0926403 36.0259904 0.30053120
    4. Comme vres n'était pas donnée, il fallait partir des valeurs expérimentales et calculer
      (tmp(2)-mres)/(tmp(3)/sqrt(numcli)) $ \approx1.32$. Ici encore, il y a un bon accord théorie-expérience.
  7. Quelle serait la durée de séjour moyenne dans une file M/M/1 ayant les mêmes valeurs de $ \lambda$ et $ \mu$? Expliquer pourquoi la valeur obtenue pour la file M/Ga/1 est inférieure.
    1. Pour une file M/M/1, le temps de séjour moyen serait

      $\displaystyle \mathrm{E}\left(séjour\right)=\frac{1}{\mu-\lambda}=\frac{1}{1/marv-1/msrv}\approx639$

    2. La valeur nettement plus faible ( $ \approx471$) obtenue avec une file M/Ga/1 ayant les mêmes débits s'explique en premier lieu par le fait que la variance des services est ici $ a\, b^{2}\approx1450$ tandis qu'elle serait $ \left(a\, b\right)^{2}$ -soit 2.4 fois plus- pour des services exponentiels de même débit.
  8. Comparer le nombre de services résiduels constatés et le nombre de clients entrés dans le système. Expliquer la proportion constatée en la reliant à $ \lambda$ et $ \mu$.
    1. Il a été constaté $ 22856$ services résiduels alors que $ numcli=25005$. Soit un ratio de $ 0.914$.
    2. La durée totale simulée est $ \approx1612343$, tandis que la durée d'occupation du serveur est (total-idle) $ \approx1472502$. Soit un ratio de $ 0.913$.
    3. Le fait que serveur et clients voient la même valeur pour le taux d'activité est lié aux arrivées exponentielles. Le fait que le taux temporel (celui vu par le serveur) soit approximativement égal à $ \rho$ est une propriété générale. 

//Corrigé du DS 23a

if MSDOS then 
  cd("c:/MapleV6/Users/Douillet/decis") ;
  libname1='decis.sce';
else
  cd("~/docs/Ensait/decis/decis_ds23a") ;
  libname1='~/docs/Ensait/decis/scilab/decis.sce';
  libname2='~/docs/Ensait/waits/scilab/waits.sce';
end//if
BASErep=pwd();
exec(libname2, -1); exec(libname1, -1)
cd (BASErep)
myfig=BASErep+'/figures/';
global gap alt lesx
global durees serveur sejours attentes autres resids services

printf('---------------------------------------------------\n');

initpick()
pr=[p1,p2,p3,p4];xdel(winsid()); 
pickc()
printf(  'gap= '); printf('%5.3f  ',gap'); 
printf('\nalt= '); printf('%5d  ',alt'); printf('\n\n');
printf('%s\n', sci2exp(gen))

chances=gap;
for j=1:4 do chances(alt(j))=chances(alt(j))+1-gap(j); end
printf('\nchances= '); printf('%5.3f  ',chances'); 
printf('\n     pr= '); printf('%5.3f  ',pr'); printf('\n');

li=[1*ones(1:3),2*ones(1:8),3*ones(1:5),4*ones(1:4)];
// tmp=tabul(li,'i'); tmp=tmp(:,2)'; tmp/sum(tmp)
histplot(bor,li)
curax=gca(); curax.children(2).children.thickness=2;

[ind,obs]=dsearch(lesx,bor);
hyp=pr*N;
printf('\nvaleur'); printf('%5d ',(1:4)')
printf('\n   hyp'); printf('%5d ',round(hyp'))
printf('\n   obs'); printf('%5d ',obs'); printf('\n');

pearson=sum((hyp-obs).^2 ./ hyp);
pearson_red=(pearson-3)/sqrt(6);

prob=cdfchi('PQ',pearson,3);
printf('pearson= %5.3f, pearson_red= %5.3f, prob= %f\n',..
  pearson, pearson_red, prob);

mu=mean(li); sigma=sqrt(mean(li.^2)- mu^2); // variance(li)
qui=[1,2,3,4]; mu=meanf(qui,pr); sigma=sqrt(meanf(qui .^2,pr)-mu^2);
// variancef(qui,pr*size(li,'*'))
moy=mean(lesx); var=variance(lesx);
moy11=mean(lesx(10:20)); var11=variance(lesx(10:20));

printf('\n   mu= %f, sigma^2=%f\n', mu, sigma^2);
printf('  moy= %f,     var=%f\n', moy, var);
printf('moy11= %f,   var11=%f\n\n', moy11, var11);

for n=[N,11] do for p=[0.975, 0.995] do
  k=cdfnor('X',0,1,p,1-p); delta=k*sigma/sqrt(n);
  printf('au risque %d%%, k= %4.2f et moy(%3d) in [%6.4f,%6.4f]\n',..
    round(200*(1-p)), k, n, mu-delta, mu+delta)
end//for
end//for

printf('\n')
n=N; for p=[0.975, 0.995] do
  k=cdft('T',n-1,p,1-p); delta=k*sqrt(var/n);
  printf('N=%3d et risque %d%%, t= %4.2f et mu in [%6.4f,%6.4f]\n', ..
    n, round(200*(1-p)), k, moy-delta, moy+delta)
end//for
n=11; for p=[0.975, 0.995] do
  k=cdft('T',n-1,p,1-p); delta=k*sqrt(var11/n);
  printf('N=%3d et risque %d%%, t= %4.2f et mu in [%6.4f,%6.4f]\n', ..
    n, round(200*(1-p)), k, moy11-delta, moy11+delta)
end//for

printf('\n---------------------------------------------------\n\n');

tmp=fun2string(petitessai,'petitessai1');
tmp=strsubst(tmp, '12345','qzzq');
tmp=strsubst(tmp, '12','16');
tmp=strsubst(tmp, '18','26');
tmp=strsubst(tmp, '15','10');
tmp=strsubst(tmp, '25','19');
tmp=strsubst(tmp, 'petitessai1','petitessai1(qzzq)');
execstr(tmp); 

tmp=fun2string(petitessai,'petitessai2');
tmp=strsubst(tmp, '12345','qzzq');
tmp=strsubst(tmp, '12','14');
tmp=strsubst(tmp, '15','9');
tmp=strsubst(tmp, '25','17');
tmp=strsubst(tmp, '18','25');
tmp=strsubst(tmp, 'petitessai2','petitessai2(qzzq)');
execstr(tmp); 

petitessai1(); printf('\n'); 

N=50000; menu_MGa1();
cd(BASErep); repertoire='1s-gam3'; cd(repertoire);
mv=drawq(); drawgam();
cd(BASErep);

tmp=msscanf(mv(5,:),'%s %f %f%f'); printf('\n%s\n',mv(5,:)); 
aaa=msrv^2/vsrv, bbb=msrv/aaa
z_moy_serv=(tmp(2)-msrv)/sqrt(vsrv/numcli)
var_s2=3*aaa*(aaa+2)*bbb^4 -vsrv^2
z_var_serv=(tmp(3)^2-vsrv)/sqrt(var_s2/numcli)

tmp=msscanf(mv(6,:),'%s %f %f%f'); printf('\n%s\n',mv(6,:)); 
mmas=(aaa+1)*bbb; vmas=mmas*bbb;
printf('mmas=%f, vmas=%f\n', mmas,vmas);
vres=vmas/3+mmas^2/12;
printf('mres=%f, vres=%f, sres=%f\n', mmas/2, vres, sqrt(vres));

shortname=['durees','autres','sejour','service','resid','','walker',''];
for fig=[1:5,7] do
curfig=scf(fig); curax=gca();
curax.title.font_size=3; curax.sub_ticks=[0,0];  
curax.margins=[0.08,0.03,0.12,0.1] ;
select fig
  case 4 then curax.data_bounds(2)=4*msrv;
  case 5 then curax.data_bounds(2)=3*msrv;
end//select
set_posfig_dim(curfig.figure_size(1),curfig.figure_size(2));
xs2eps(fig, myfig+shortname(fig));
end
previous up next_inactive
Up: Return to previous menu


douillet@ensait.fr
2008-01-27