//////////////////////////////////////////////////
//Codes scilab présents dans le fichier Survie.pdf
//////////////////////////////////////////////////

//1 Générer des variables aléatoires

//1.1 La fonction grand
    X=grand(23,1,'def');
    help grand;

//1.2 La simulation de lois plus générales
    //à la main
    P=[1/6; 1/4; 1/6; 1/6; 1/6; 1/12];
    I=cumsum(P);//on définit les bornes sup de chaque sous-intervalle
    U=grand(1,1,'def');
    X=min(find(U<I));//find donne les indices de I qui vérifient l'inégalité  
    
    //utilisation de la loi multinomiale
	P=[1/6; 1/4; 1/6; 1/6; 1/6];//le vecteur est de dimension 5!
	grand(2,'mul',1,P);//deux vecteurs colonnes élémentaires
	//le 1 donnant la face sortie à chaque tour

//1.3 Les processus classiques
    
    //Chaîne de Markov homogène
    P=[0.1 0.3 0 0.6;
	0.4 0.4 0.1 0.1;
	0 0.8 0 0.2;
	0.2 0 0.2 0.6];
	x0=1;
	X=grand(10,'markov',P,x0);
    
    //définition d'une matrice diagonale
    D=[1:10];
    D=diag(D);
    
    //génération d'une trajectoire de la ruine du joueur
    function [X]=generateur(f,b,p,nb_pas);//0<f<b entiers 
		z=zeros(1,b+1);
		u=z;u(2)=1-p;
		v=z;v(2)=p;
        P=toeplitz(u,v);
        P(1,1)=1;P(b+1,b+1)=1;//états absorbants
        P(1,2)=0;P(b+1,b)=0;//on retire les valeurs inutiles
        disp(P);
		X=grand(nb_pas,'markov',P,f+1);
		if X($)==b+1 then
			disp('le joueur remporte le jackpot');
		elseif X($)==1
			disp('le joueur est ruiné');
        else disp('il reste au joueur la somme de '+string(X($)-1));
		end
	endfunction
    
    //Chaîne plus générale
    function [X]=generateur(N,x0);//N=nb de pas, x0=pos. initiale
		X=[x0];//initialisation de la marche
		last_position=x0;//derniere visite de la chaîne
		for k=1:N
			last_position=grand(1,1,'nor',last_position,1);
			X=[X last_position];
		end;
	endfunction;		
    
    //Processus de Poisson
    //génération d'un processus de Poisson(lambda) sur [0,t]
	function [T]=generateur(t,lambda);
		while T($)<t
			E=grand(1,1,'exp',1/lambda)
			T=[T T($)+E];
		end;
	endfunction;

//2 Illustrations

//2.1 Faire un graphique
    
    //La fonction constante par morceaux
    x=[0:7];
	P=[1/6; 1/4; 1/6; 1/6; 1/6; 1/12];
	P=cumsum(P);//vecteur des sommes partielles
	P=[0; P; P($)];
	clf;//efface les tracés précédents
	plot2d2(x,P,5);

    //La fonction lissée
    x=[-7:.01:7];
	y=x-3;
	y=exp(-y^{2}/4);
	y=1/(2*sqrt(2*%pi))*y;
	clf;
	plot2d(x,y,2);

    //L'histogramme
    X=grand(10000,1,'exp',.5);
	clf;
	histplot(20,X,9);
	x=[0:.01:max(X)];
	plot2d(x,2*exp(-2*x),5);

    //Les tracés 2d de manière générale
    function [X]=trajectoire(nb_pas,P);//P proba sur 1,2,3,4
        X=[0; 0 ];
        nord=[0; 1];
        sud=[0; -1];
        est=[1; 0];
        ouest=[-1; 0];
        M=[nord sud est ouest];
        P(4)=[];//se souvenir que l'on ne doit avoir que 3 coordonnées
        for k=1:nb_pas
            sens=grand(1,'mul',1,P);
            X(:,k+1)=X(:,k)+M*sens;            
        end
        plot2d(X(1,:),X(2,:),2);
    endfunction

    //La décoration
    X=grand(10000,1,'exp',.5);
	clf;
	histplot(20,X,9);
	x=[0:.01:max(X)];
	plot2d(x,2*exp(-2*x),5);
	title('Illustration de la convergence vers une densité','fontsize',4);
	xlabel('Axe des x','fontsize',4);
	ylabel('Axe des y','fontsize',4);
	legends(['Histogramme';'Densité'],[9;5],opt='ur', font_size=4);

    //La multiplication des graphes dans une même fenêtre
    X=grand(10000,1,'exp',.5);
	clf;
    subplot(1,2,1)//premier graphe
	histplot(20,X,9);
	x=[0:.01:max(X)];
	plot2d(x,2*exp(-2*x),5);
	title('Illustration de la convergence vers une densité','fontsize',2);
	xlabel('Axe des x','fontsize',2);
	ylabel('Axe des y','fontsize',2);
	legends(['Histogramme';'Densité'],[9;5],opt='ur', font_size=2);
    subplot(1,2,2)//deuxième graphe
    x=[-7:.01:7];
	y=x-3;
	y=exp(-y^{2}/4);
	y=1/(2*sqrt(2*%pi))*y;
	plot2d(x,y,2);    
	title('Densité de la gaussienne de moyenne 3 et de variance 4','fontsize',2);
	xlabel('Axe des x','fontsize',2);
	ylabel('Axe des y','fontsize',2);
    
//2.2 Ce qu'il faut y mettre

    //La convergence en loi
    function illustration();
		x=grand(1,100,'geom',0.2);
		y=grand(1,1000,'geom',0.2);
		z=grand(1,10000,'geom',0.2);
		x=tabul(x,'i');//i pour increasing
		y=tabul(y,'i');//le tri par défaut est décroissant!
		z=tabul(z,'i');
		x(:,2)=x(:,2)/100;//normalisation des occurences
		y(:,2)=y(:,2)/1000;
		z(:,2)=z(:,2)/10000;
        m=max(x($,1),y($,1),z($,1));
		clf;
		//on réprésente alors les fonctions de répartition
		//empiriques, avec la fonction cumsum
        xx=[0; x(:,1) ;m];xxx=[0; x(:,2) ; x($,2)];
        yy=[0; y(:,1) ;m];yyy=[0; y(:,2) ; y($,2)];
        zz=[0; z(:,1) ;m];zzz=[0; z(:,2) ; z($,2)];
		plot2d2(xx,cumsum(xxx),3);
		plot2d2(yy,cumsum(yyy),4);
		plot2d2(zz,cumsum(zzz),5);
		title('Illustration de la convergence des fonctions de répartition empiriques','fontsize',4);
		xlabel('Axe des x','fontsize',4);
		ylabel('Axe des y','fontsize',4);
		legends(['n=100';'n=1000';'n=10000'],[3;4;5],opt='ur', font_size=4);
	endfunction;

    //L'intervalle de confiance
    function tcl(nb_iterations);
        X=grand(1,nb_iterations,'exp',5);
        n=[1:nb_iterations];
        X=cumsum(X)./n;
        ecart=5*1.96./sqrt(n);
        clf;
        plot2d2(n',[X'-ecart' X' X'+ecart'],[5; 2; 5]);
        plot2d(n,ones(n)*5);
		title('Illustration du TCL','fontsize',4);
		xlabel('Taille de l''échantillon','fontsize',4);
		ylabel('Axe des y','fontsize',4);
		legends(['borne haute de l''IC';'estimation de la moyenne';'borne basse de l''IC';'moyenne théorique'],[5;2;5;1],opt='ur', font_size=4);
    endfunction

//3 Tests statistiques

//3.2 Ce que l'on peut faire faire à Scilab lors d'un test

    //A la recherche des quantiles
    P=cdfnor('PQ',5,2,2);//évaluation de la fonction de répartition au point X=5
    //pour une gaussienne (2,4)
    X=cdfnor('X',0,1,.975,.025);//quantile associé à P=.975
    //pour une normale centrée réduite
    
    //L'apprentissage par l'exemple
    s=cdfnor('X',75,sqrt(.1),.99,.01);
    x_obs=[75.2 72.6 74.5 75 75.5 73.7 74.1 75.8 75.8 75];
    T_obs=mean(x_obs);
    
    //test dans le cas général
    function test(x_obs);
		n=length(x_obs);
		T_obs=mean(x_obs);
		s=cdfnor('X',75,sqrt(1/n),.01,.99);
		decision=1*(T_obs<= s);
		disp('H'+string(decision));
	endfunction;

//4 Identifier une erreur

function generateur();
    X=[0];
    for k=1:10
        disp(k);
        X(k)=1/(10-k);
    end
endfunction