//-------------------------------------- // Université de Nice -- Sophia Antipolis // Polytech'Nice -- Sophia // Département Électronique, 3e année // Cours de Statistique Appliquée // Iannis.Aliferis@unice.fr //--------------------------------------- // Ecécuter le script à partir du fichier // > exec('le_nom_du_fichier.sce' ) // // Vous pouvez, à la fin, obtenir les valeurs // de toutes les variables calculées en tapant leur nom // dans la ligne de commande. //--------------------------------------- // Chapitre : // Statistique Inférentielle: échantillonnage // TD 6, Théorie d'échantillonnage, simulations avec Scilab //--------------------------------------- // Script : // Explorer le comportement de la moyenne d'un échantillon; // comparer l'expérience à la théorie //--------------------------------------- // Fonction qui change le titre d'une figure de la // maniere suivante : // fenetre : numero de la fenetre sur laquelle on travaille // xlim : limites des abcisses sous la forme [x_inf,x_sup] // ylim : idem pour les ordonnées // si on ne veut pas changer xlim et ylim, passer [] comme argument // labelx : le texte associe aux abscisses // labely : le texte associe aux ordonnees function titleit(fenetre,xlim,ylim,labelx,labely,titre) scf(fenetre); axes=get("current_axes"); if (xlim ~= []) axes.data_bounds=[xlim,ylim]; end axes.x_label.text=labelx; axes.y_label.text=labely; axes.title.text=titre; endfunction //-------------------------------------- // Générer la population //-------------------------------------- // On définit ici une population de taille finie, égale à N N = 1e6; // Si on veut une population de taille infinie, on utilise directement cdfnor() etc // dans la boucle d'échantillonnage (cf. commandes en commenaire dans la boucle) // La loi p_X(x) suivie par la population: normale, uniforme ou autre. // Penser ne pas regarder que la loi normale! // population normale mu_X = 0 sigma_X = 1 population = rand( N, 1,'normal')*sigma_X+mu_X; // A VOUS : si le temps le permet, essayez avec d'autres populations // par exemple uniforme // population uniforme //a = 0; //b = 1; //mu_X = (a+b)/2; //sigma_X = sqrt( (a-b)^2 / 12 ); //population = a+rand( N,1,'uniform')*(b-a); //-------------------------------------- // Créer les vecteurs x.bar.vec et z.exp.vec //-------------------------------------- // nombre d'échantillons à prélever k = 10e3; x_bar_vec = zeros(k,1); z_exp_vec = zeros(k,1); //-------------------------------------- // Prélever k échantillons, de taille n, de la population //-------------------------------------- // 6.1.4 : étude paramétrique : on va créer une boucle sur différentes tailles d'échantillon // taille de l'échantillon n = [3,10,30,50]; // variables qui contiennent les grandeurs desirees norm_ecart_rel_freq.rel = zeros(length(n),1); max_abs_ecart_rel_freq_rel = zeros(length(n) ,1); which_max_abs_ecart_rel_freq_rel = zeros( length(n),1 ); for (index_n=1:length(n)) for (i =1:k) // avec ou sans remplacement, ça n'a pas trop d'importance si n << N index=N*rand(n(index_n),1,"uniform"); echantillon=population(ceil(index)); x_bar_vec(i) = mean( echantillon ); z_exp_vec(i) = ( x_bar_vec(i) - mu_X ) / ( sigma_X / sqrt(n(index_n)) ); end //-------------------------------------- // Afficher l'histogramme de Z //-------------------------------------- // on montre la densité, c-à-d la fréquence relative dans chaque classe est divisée par la largeur de la classe. scf(0); clf; histplot( 11,z_exp_vec); // superposer la ddp théorique de Z: z_th_vec = min(z_exp_vec):(max(z_exp_vec)-min(z_exp_vec))/500:max(z_exp_vec) ; pz_theorique=1/sqrt(2*%pi)*exp(-(z_th_vec).^2/2); plot( z_th_vec, pz_theorique); titre=strcat(['Densite de probabilite, n=',string(n(index_n))]); titleit(0,[],[],'z_exp','densite',titre); //---------------------------------------------- // 6.1.2 Approche quantitative //---------------------------------------------- // partager l'histogramme en M+1 classes M = 99; // calculer les M quantiles de la ddp théorique // les probabilités associées p = 1/(M+1):1/(M+1) : 1 - 1/(M+1); // et les quantiles correspondant à ces probabilités // A VOUS : utilisez cdfnor // qu'on utilise pour fixer les bornes (breaks) des classes bornes = [ min(z_exp_vec), q_m_th, max(z_exp_vec) ]; // tracer l'histogramme avec les classes précises scf(1); clf; histplot(bornes, z_exp_vec); // et superposer la ddp théorique plot( z_th_vec, pz_theorique); titre=strcat(['Densite de probabilite, n=',string(n(index_n))]); titleit(1,[],[],'z_exp_vec','densite',titre); // tracer l'histogramme avec les classes précises, *en fréquences* // cela devrait donner la même fréquence pour toutes les classes, // puisque on a défini des classes "équiprobables": // les k valeurs devraient être équiréparties dans les M+1 classes, // ce qui donne une fréquence théorique égale à k/(M+1) // (R donne des warnings parce que on utilise des fréquences avec des largeurs inégales!) scf(2); clf; // A VOUS : histplot avec les bons arguments // ajouter une ligne horizontale au niveau de la // fréquence théorique que toutes les classes //devraient présenter // A VOUS : plot avec les bons arguments titre=strcat(['frequences non normalisees, n=',string(n(index_n))]); titleit(2,[],[],'z_exp_vec','frequences',titre); // les fréquences relatives expérimentales (à partir de l'histogramme) for freqi=1:length(bornes)-1 freq_rel_exp(freqi)=sum((z_exp_vec>bornes(freqi))&(z_exp_vec