//--------------------------------------
// 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<bornes(freqi+1)))/k;
	end

	// les fréquences relatives théoriques
	// on a construit M+1 classes "équiprobables"
	freq_rel_th = 1 / ( M+1 );

	// l'écart relatif entre les fréquences relatives expérimentales et théoriques
	ecart_rel_freq_rel = ( freq_rel_exp - freq_rel_th ) / freq_rel_th;
      
	// afficher les valeurs absolues de l'écart relatif
	scf(3);
	clf;
	plot( abs( ecart_rel_freq_rel ));
        titre=strcat(['valeurs absolues de l ecart relatif, n=',string(n(index_n))]);
        titleit(3,[0,M+1],[0,1],'index','ecart_rel_freq_rel',titre);
 
	norm_ecart_rel_freq_rel(index_n) = sqrt( sum( ecart_rel_freq_rel^2 ) );
        [mmm,www]=max(abs(ecart_rel_freq_rel));     
	max_abs_ecart_rel_freq_rel(index_n) = mmm;
	which_max_abs_ecart_rel_freq_rel(index_n) = www;

//----------------------------------------------
// 5.1.3	  Quantiles vs. Quantiles
//----------------------------------------------
      // quantiles théoriques : q_m_th
      // quantiles expérimentaux : q_m_exp

      // A VOUS : trouvez les quantiles experimentaux via  la commande
      // perctl()
      scf(4);
      clf;
      plot(q_m_th,q_m_exp,'o');
      // A VOUS : trouvez la droite de regression entre q_m_th et q_m_exp 
      // reglin()
      xx=[-2.5:2.5];
      reg=a*xx+b;
      plot(xx,reg);

	titre=strcat(['Quantiles theoriques vs quantiles exper., n=',string(n(index_n))]);
      titleit(4,[],[],'quantiles theoriques','quantiles experimentaux',titre);
 

	// tracer la droite théorique
	// A VOUS ( via plot())

pause;
end// end loop on n, sample size

//----------------------------------------------
// 6.1.4 : diagrammes de l'étude paramétrique
//----------------------------------------------
scf(5);
clf;

plot( n, norm_ecart_rel_freq_rel )
scf(6);
clf;
plot( n, max_abs_ecart_rel_freq_rel )
scf(7);
clf;
plot( n, which_max_abs_ecart_rel_freq_rel )



