Matlab

De Philippedia
Aller à : navigation, rechercher

Pour Philippe, Matlab est un langage sans émotion. Qui, paradoxalement, peut donner d'intenses frustrations, ou de grandes joies de soulagement ("ça marche enfin!"). Finalement, Philippe n'a plus vraiment eu l'occasion de s'y plonger. Le projet 2 de risques hydrologiques est passé! (merci Steph!)

D'ailleurs, le voici ! :) (merci Camille, Jonas et Antoine!)


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Chargement des données
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  data = load('StatData/dataP.txt');
  pluie = data(:,4)/10;%en mm
  an=data(:,1);
  %hist(pluie,100)
  %ecdf(pluie)
  %mean(pluie)
  %var(pluie)

  %plot(pluie,'.')
  pause on


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Problème 1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Question a)

  % Détermination des cumuls annuels de pluie (en mm)
  for i = 1:25
     annee = 1973 + i ;
     label_annee(i)=annee;
     liste = find(data(:,1)==annee) ;
     pluieAnnee(i)  = sum(data(liste,4)/10) ;
     %% find(ARG) donne les indexes pour lesquelle ARG est valable
  end;

  bar(label_annee,pluieAnnee);
  hold on
  line([1974,1998],[mean(pluieAnnee),mean(pluieAnnee)]);
  moy_annuel=mean(pluieAnnee)
  sigma_an=sqrt(var(pluieAnnee))

  pause
  %Détermination de la distribution des cumuls mensuels

  for i = 1:25
     annee = 1973 + i ;
     pliste = find(data(:,1)==annee) ;
     for j = 1:12
        liste = find(data(pliste,2)==j) ;
        pluieMois((i-1)*12+j,1) = annee;
        pluieMois((i-1)*12+j,2) = j;
        pluieMois((i-1)*12+j,3) = sum(pluie(pliste(liste))) ;
        %% find(ARG) donne les indexes pour lesquelle ARG est valable
     end;
  end;

  hold off
  hist(pluieMois(:,3),0:10:310) %histogramme
  [P1, x_pluie] = ecdf(pluieMois(:,3));
  m_pluie = mean(pluieMois(:,3))
  sigma=sqrt(var(pluieMois(:,3)))

  pause
  hold off
  ecdf(pluieMois(:,3)) % courbe de distribution cumulée

  %pause
  %hold on
  %plot(cdf('exp',x_pluie,1/m_pluie),P1)

  pause

  %Détermination des moyennes mensuelles
  for i = 1:12
    liste = find(pluieMois(:,2)==i) ;
    pluieMoisMoy(i) = mean(pluieMois(liste,3)) ;
    pluieMoisVar(i) = sqrt(var(pluieMois(liste,3)));
    %% find(ARG) donne les indexes pour lesquelle ARG est valable
 end;

 hold off
 bar(pluieMoisMoy);
 hold on
 line([1,12],[mean(pluieMoisMoy),mean(pluieMoisMoy)]);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Question b)

 %Pluies journalière dépassant 1mm
 j=1;
 for i=1:length(pluie)
     if pluie(i)>=1
         pluie1mm(j)=pluie(i);
         j=j+1;
     end
 end

 %Densité de probabilité empirique
 hold off
 hist(pluie1mm,0:1:120)
 pause
 %Fonction de répartition empirique
 ecdf(pluie1mm);

 %Ajustement exponentielle
 moy=mean(pluie1mm)
 landa=1/moy

 %intevalle de confiance sur landa
 n=length(pluie1mm);
 z=1.95996;
 I=[landa-z*1/sqrt(n) landa+z*1/sqrt(n)]

 %a=1
 %[int] = grpstats(pluie1mm,a,'meanci')

 %[m,p,g] = grpstats(Weight,Model_Year,{'mean','predci','gname'})
 %n = length(m)
 %errorbar((1:n)',m,p(:,2)-m)
 %set(gca,'xtick',1:n,'xticklabel',g)
 %title('95% prediction intervals for mean weight by year')

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Question c)

 % voir avec toolbox :statistics: distribution fitting tool

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Question d) stationarité

 plot(1:length(pluie),cumsum(pluie>=1))
 hold on
 plot([1,length(pluie)],[1,sum(pluie>=1)],'r')
 hold off
 %equation de la droite: y = p1*x + p2, trouvée avec figure, tool, basic
 %fitting
   p1 = 0.33866;
   p2 = 0.66134;
 %calcul des résidus
 xr=1:length(pluie);
 yr=cumsum(pluie>=1)-(p1*xr + p2)';
 mean(yr)
 plot(xr,yr)
 hold on
 plot([1,length(pluie)],[0,0])

 pause
 plot(an,cumsum(pluie>=1))
 hold on
 plot([min(an),max(an)],[1,sum(pluie>=1)],'r')
 hold off

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Question e) nombre de chute de pluie dépassant un seuil par année

 %Pluies journalière dépassant un seuil s
 s=1 %mm

 j=1;
 for i = 1:25
    annee = 1973 + i ;
    pliste = find(data(:,1)==annee) ;

 for y=1:length(pliste)
     if pluie(pliste(y))>=s
         pluieseuil1(j,1)=annee;
         pluieseuil1(j,2)=pluie(pliste(y));
         j=j+1;
     end
 end
 end

 pluieseuil1;

 j=1;
 for i = 1:25
    annee = 1973 + i ;
    label_annee(i)=annee;
    pliste = find(pluieseuil1(:,1)==annee) ;
    nbrsup1(i)=length(pliste);
 end

 bar(label_annee,nbrsup1)

 pause
 s=20 %mm

 j=1;
 for i = 1:25
    annee = 1973 + i ;
    pliste = find(data(:,1)==annee) ;

 for y=1:length(pliste)
     if pluie(pliste(y))>=s
         pluieseuil20(j,1)=annee;
         pluieseuil20(j,2)=pluie(pliste(y));
         j=j+1;
     end
 end
 end

 pluieseuil20;

 j=1;
 for i = 1:25
    annee = 1973 + i ;
    label_annee(i)=annee;
    pliste = find(pluieseuil20(:,1)==annee) ;
    nbrsup20(i)=length(pliste);
 end

 bar(label_annee,nbrsup20)

 hist(nbrsup1);
 pause
 hist(nbrsup20);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %f) Ajustement loi de poisson: paramètre landa et intervalle de confiance à
 %95% avec le maximum de vraissemblance

 [landa1,int951] = poissfit(nbrsup1)
 [landa20,int9520] = poissfit(nbrsup20)
 %méthode des moments: landa=moyenne /plot avec le tool box
 mean(nbrsup1)
 mean(nbrsup20)

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %g)log vraisemblance pour s=20mm
 landa=landa20
 n=length(nbrsup20)
 somme1=0;
 somme2=0;
 for i=1:n
     somme1=somme1+log(factorial(nbrsup20(i)));
     somme2=somme2+nbrsup20(i);
 end

 l=log(landa)*somme2-somme1-landa*n

 landaplot=12.42:0.001:12.62;
 lplot=log(landaplot)*somme2-somme1-landaplot*n;
 plot(landaplot, lplot)

 j = find(lplot == max(lplot));
 landaplot(j)

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %h) méthode Bayesienne
 x=nbrsup20
 %f(x|lambda):
 f=1;
 for i=1:25
 p=poisspdf(x(i),landa);
 f=f*p;
 end
 f
 %pi(lambda)
 Y = unifpdf(landa,5,20)
 %Q=f*pi
 Q=f*Y


 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %Problème 2
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%
 % A
 %%%%%%%%%
 clear all ;

 data = load('StatData/dataP.txt');
 pluie = data(:,4)/10;%en mm
 an=data(:,1);

 debut = min(an) ; % prmière année
 fin   = max(an) ; % deriniere annee
 j=1 ;
 for i = debut : fin
     indexes             = find(an == i) ;
     maxima_annuels(j)   = max(pluie(indexes)) ;
     j = j + 1 ;
 end;

 bar(min(an):max(an),maxima_annuels)
 title('Maxima annuel des précipitations')
 ylabel('Pluie [mm]')
 xlabel('Année')
 pause

 hist(maxima_annuels,5)
 ylabel('Nombre dannées')
 xlabel('Précipitation maximale annuelle [mm]')
 pause

 [parametres, parametresIC] = gevfit(maxima_annuels) ;

 parametres    % Parametre xi, sigma, mu
 parametresIC   % intervalles de confiance à 95%
 pause

 dataSort = sort(maxima_annuels);
 for i = 1:length(dataSort)
     F(i) = (i-0.5)/(length(dataSort));
 end

 stairs(dataSort,F), hold on
 plot(sort(maxima_annuels), gevcdf(sort(maxima_annuels),parametres(1),parametres(2),parametres(3)),'g')
 plot(sort(maxima_annuels), gevcdf(sort(maxima_annuels),parametresIC(1,1),parametresIC(1,2),parametresIC(1,3)),'r')
 plot(sort(maxima_annuels), gevcdf(sort(maxima_annuels),parametresIC(2,1),parametresIC(2,2),parametresIC(2,3)),'r'), hold off
 title('Fonction de probabilité empirique et estimée')
 xlabel('Précipitation maximale annuelle [mm]')
 legend('probabilité empirique','probabilité estimée','intervalle de confiance 95%',4)
 pause

 [parametresGumbel, parametresGumbelIC] = evfit(maxima_annuels) ;

 parametresGumbel    % Parametre mu sigma
 parametresGumbelIC   % intervalles de confiance à 95%
 pause

 plot(sort(maxima_annuels),sort(maxima_annuels)), hold on
 plot(gevinv(F,parametres(1),parametres(2),parametres(3)),sort(maxima_annuels),'.'), hold off
 xlabel('Quantile selon Fréchet [mm]')
 ylabel('Quantile empirique [mm]')
 pause

 %%%%%%%%%
 % B
 %%%%%%%%%

 tps_r = 2:1000;
 tps_r_empiriques = 1./(1-F);
 quantile_C = gevinv((1-1./tps_r),parametres(1),parametres(2),parametres(3));
 quantile_Cplus = gevinv((1-1./tps_r),parametresIC(2,1),parametresIC(2,2),parametresIC(2,3));
 quantile_Cmoins = gevinv((1-1./tps_r),parametresIC(1,1),parametresIC(1,2),parametresIC(1,3));

 semilogx(tps_r_empiriques,dataSort,'.'), hold on
 semilogx(tps_r,quantile_C)
 plot([100 100],[0 300],'k')
 plot([1 100],[quantile_C(99) quantile_C(99)],'k:'), hold off
 axis([1 1000 0 300])
 legend('Quantiles empiriques','Quantiles selon Féchet',2)
 xlabel('Période de retour [an]')
 ylabel('Pluies [mm]')
 pause

 semilogx(tps_r_empiriques,dataSort,'.'), hold on
 semilogx(tps_r,quantile_C)
 semilogx(tps_r,quantile_Cplus,'r')
 semilogx(tps_r,quantile_Cmoins,'r')
 semilogx(tps_r_empiriques,dataSort,'.')
 plot([100 100],[0 800],'k')
 plot([1 100],[quantile_C(99) quantile_C(99)],'k:')
 plot([1 100],[quantile_Cplus(99) quantile_Cplus(99)],'k:')
 plot([1 100],[quantile_Cmoins(99) quantile_Cmoins(99)],'k:'), hold off
 axis([1 1000 0 800])
 legend('Quantiles empiriques','Quantiles selon Féchet','Quantiles "limites" selon Féchet 95%',2)
 xlabel('Période de retour [an]')
 ylabel('Pluies [mm]')
 pause

 [ quantile_C(99) quantile_Cplus(99) quantile_Cmoins(99)]
 pause

 %%%%%%%%%
 % C
 %%%%%%%%%

 f1 = figure ; hold on ;
 f2 = figure ; hold on ;

 %% Espérance conditionnelle
 for i = 1:100
     s(i)=max(pluie)/100 * i;
     Y = pluie(find(pluie>s(i))) ;
     E_s(i) = mean(Y-s(i)) ;
 end

 % Choix "visuel"
 clear X;
 clear Y;
 X(1) = 0 ; X(2) = 60;
 Y(1) = 8 ; Y(2) = 22;

 plot(s,E_s,'.k'), hold on
 plot(X,Y,'r'), hold off

 [((Y(2)-Y(1))/(X(2)-X(1))) Y(1)] %xi, sigma

 figure(f1) ; plot(s,E_s,'.k')

 %% Définition du seuil

 choix_seuil = [10 20 30 40]

 for k = 1 : length(choix_seuil)

     seuil = choix_seuil(k)

     clear Y
     Y   = pluie(find(pluie>seuil)) -seuil ;
     n_s = length(Y) ;
     n_a = length(pluie) ;

     [gpparam , gpparam_ci] = gpfit(Y)

     xi_pareto = gpparam(1) ; sigma_pareto = gpparam(2) ;

     x(1) = seuil ; x(2) = max(s) ;
     y = (sigma_pareto + xi_pareto * (x-seuil) )/(1-xi_pareto) ;

     figure(f1)
     plot( x , y, '-r') ;
     text( x(2), y(2) , sprintf('seuil %3.1f mm \\rightarrow',seuil),'HorizontalAlignment','right')  ;
      xlabel('seuil') ;
      ylabel('E[X-s|X>s]') ;

     figure(f2)
     lambda = length(pluie)/length(Y) / 365;
     plot(     sort(Y+seuil)   ,  1./(1-gpcdf(sort(Y),xi_pareto,sigma_pareto))*lambda       ,'-r') ;
     text( max(sort(Y+seuil))  ,  1./(1-gpcdf(max(sort(Y)),xi_pareto,sigma_pareto))*lambda  ,...
          sprintf('seuil %3.1f mm \\rightarrow',seuil),'HorizontalAlignment','right')
      xlabel('Précipitation en mm') ;
      ylabel('Temps de retour en années') ;
 end
 pause

 % pour la suite nous choisissons un seuil de 30 mm, bien représentatif du
 % "faisseau" 10-40.

     seuil = 30;

     clear Y
     Y   = pluie(find(pluie>seuil)) -seuil ;
     n_s = length(Y) ;
     n_a = length(pluie) ;

     [gpparam , gpparam_ci] = gpfit(Y)

     xi_pareto = gpparam(1) ; sigma_pareto = gpparam(2) ;

     x(1) = seuil ; x(2) = max(s) ;
     y = (sigma_pareto + xi_pareto * (x-seuil) )/(1-xi_pareto) ;

     plot(s,E_s,'.k'), hold on
     plot( x , y, '-r'), hold off
      xlabel('seuil') ;
      ylabel('E[X-s|X>s]') ;
     pause

     lambda = length(pluie)/length(Y) / 365;

     plot(sort(Y+seuil),1./(1-gpcdf(sort(Y),xi_pareto,sigma_pareto))*lambda,'-r')
       xlabel('Précipitation en mm') ;
      ylabel('Temps de retour en années') ;
      pause

 dataSort = sort(Y);
 for i = 1:length(dataSort)
     F(i) = (i-0.5)/(length(dataSort));
 end

      plot(dataSort+seuil,dataSort+seuil), hold on
 plot(gpinv(F,xi_pareto,sigma_pareto)+seuil,dataSort+seuil ,'.'), hold off
 xlabel('Quantile selon Pareto [mm]')
 ylabel('Quantile empirique [mm]')
 pause

 % pour la pseudo-période de retour, nous avons :

 m_r = 2:1000;
 m_r_empiriques = 1./(1-F);
 quantile_C_pareto = gpinv((1-1./m_r),xi_pareto,sigma_pareto)+seuil;
 pause

 semilogx(m_r_empiriques,dataSort+seuil,'.'), hold on
 semilogx(m_r,quantile_C_pareto), hold off
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Pseudo-période de retour m')
 ylabel('Pluies [mm]')
 pause

 % On converti cette loi en loi des valeurs extrêmes

 xi = xi_pareto;
 mu = seuil + sigma_pareto/xi_pareto*((n_s/25)^xi_pareto - 1);
 sigma = sigma_pareto * (n_s/25)^xi_pareto;

 xi_plus = gpparam_ci(2,1);
 mu_plus = seuil + gpparam_ci(2,2)/gpparam_ci(2,1)*((n_s/25)^gpparam_ci(2,1) - 1);
 sigma_plus = gpparam_ci(2,2) * (n_s/25)^gpparam_ci(2,1);

 xi_moins = gpparam_ci(1,1);
 mu_moins = seuil + gpparam_ci(1,2)/gpparam_ci(1,1)*((n_s/25)^gpparam_ci(1,1) - 1);
 sigma_moins = gpparam_ci(1,2) * (n_s/25)^gpparam_ci(1,1);

 [xi sigma mu]
 [xi_plus sigma_plus mu_plus]
 [xi_moins sigma_moins mu_moins]
 pause

 tps_r = 2:1000;
 quantile_C_pareto = gevinv((1-1./tps_r),xi,sigma,mu);
 quantile_Cplus_pareto = gevinv((1-1./tps_r),xi_plus,sigma_plus,mu_plus);
 quantile_Cmoins_pareto = gevinv((1-1./tps_r),xi_moins,sigma_moins,mu_moins);

 semilogx(m_r_empiriques*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_r,quantile_C_pareto)
 plot([100 100],[0 300],'k')
 plot([0.1 100],[quantile_C_pareto(99) quantile_C_pareto(99)],'k:'), hold off
 axis([0.2 1000 0 300])
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Période de retour [an]')
 ylabel('Pluies [mm]')
 pause

 semilogx(m_r_empiriques*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_r,quantile_C_pareto)
 semilogx(tps_r,quantile_Cplus_pareto,'r')
 semilogx(tps_r,quantile_Cmoins_pareto,'r')
 plot([100 100],[0 800],'k')
 plot([0.1 100],[quantile_C_pareto(99) quantile_C_pareto(99)],'k:')
 plot([0.1 100],[quantile_Cplus_pareto(99) quantile_Cplus_pareto(99)],'k:')
 plot([0.1 100],[quantile_Cmoins_pareto(99) quantile_Cmoins_pareto(99)],'k:'), hold off
 axis([0.2 1000 0 800])
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Période de retour [an]')
 ylabel('Pluies [mm]')
 pause

 [ quantile_C_pareto(99) quantile_Cplus_pareto(99) quantile_Cmoins_pareto(99)]
 pause

 %Comparaison Frechet - Pareto :

 semilogx(m_r_empiriques*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_r,quantile_C_pareto)
 semilogx(tps_r,quantile_C,'r')
 plot([100 100],[0 300],'k')
 plot([0.1 100],[quantile_C(99) quantile_C(99)],'k:')
 plot([0.1 100],[quantile_C_pareto(99) quantile_C_pareto(99)],'k:'), hold off
 axis([0.2 1000 0 300])
 legend('Quantiles empiriques','Quantiles selon Pareto','Quantiles selon Fréchet',2)
 xlabel('Période de retour [an]')
 ylabel('Pluies [mm]')


 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %Problème 3
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%
 % A
 %%%%%%%%%

 clear F
 clear u
 for i = 1:length(sort(maxima_annuels))
     F(i) = (i-0.5)/(length(maxima_annuels));
     u(i) = -log(-log(F(i)));
 end

 plot(u,sort(maxima_annuels),'.'), hold on
 paramG = polyfit(u,sort(maxima_annuels),1);
 Gp = paramG(1);
 ordonnee = paramG(2);
 clear Y
 Y(1) = min(u)*Gp + ordonnee;
 Y(2) = max(u)*Gp + ordonnee;
 plot([min(u) max(u)],Y,'r'), hold off
 xlabel('u = -ln(-ln(F))')
 ylabel('Pluies maximales annuelles [mm]')
 text(-1,120,sprintf('pente %1.3f',Gp))

 S = 77.8;     % km2
 i = 8;        % pente
 L = 22000;    % m
 R = 1    ;    % coefficient régional

 T2 = 100;     % ans
 T1 = 10;      % ans
 P2 = 167;     % mm, calculé avec Pareto
 P1 = 97.5;    % mm, avec Pareto

 Q10 = 28.7;   % m3/s selon Weibull (dans la théorie des valeurs extrêmes)

 Q10c = R*S^(0.8)*(P1/80)^2; %m3/s avec la méthode du crupédix

 % Avec la formule de Turraza, on peut connaître le temps de concentration

 tc1 = 0.108*((S*L)^(1/3))/i^(1/2);% en heures

 % Avec la formule de Kirpich, on peut connaître le temps de concentration

 tc2 = 19.47 * 10^(-3)*L^(0.77)/i^(0.385)/60;

 % Avec la formule de Kirpich, on peut connaître le temps de concentration

 tc3 = 145*S^(0.312)/i^(0.325)/60;

 tc = 5 %heures

 % Pour des données suffisamment longue, on peut considérer que le gradex
 % des débits équivaut à celui des pluies exprimés dans les mêmes unités

 Gq = (S/(3.6*tc))*Gp;       %le gradex des débits

 Q100 = Q10 + Gq*log(T2/10);     %débit centennal calculé avec Q10 Weibull

 Q100c = Q10c + Gq*log(T2/10);   %débit centennal calculé avec Q10 crupédix

 Qs10 = Q10c/S

 %%%%%%%%%
 % B
 %%%%%%%%%
 dataQ = load('StatData/dataQ.txt');
 debit = dataQ(:,4); %en m3/s
 an=dataQ(:,1);

 % moyenne des cumul des débit mensuels
 for i = 1:25
    annee = 1973 + i ;
    pliste = find(dataQ(:,1)==annee) ;
    for j = 1:12
        liste = find(dataQ(pliste,2)==j) ;
        debitMois((i-1)*12+j,1) = annee;
        debitMois((i-1)*12+j,2) = j;
        debitMois((i-1)*12+j,3) = sum(debit(pliste(liste))) ;
        %% find(ARG) donne les indexes pour lesquelle ARG est valable
    end;
 end;

 for j = 1:12
        debitMoisMoy(j) = mean(debitMois(12*(0:24)+j,3)) ;
 end;
 bar(debitMoisMoy)

 % Maximas annuels

 debut = min(an) ; % prmière année
 fin   = max(an) ; % deriniere annee
 j=1 ;
 for i = debut : fin
     indexes             = find(an == i) ;
     maximaQ_annuels(j)   = max(debit(indexes)) ;
     j = j + 1 ;
 end;

 pause

 %Détermination de la distribution des cumuls mensuels

 bar(min(an):max(an),maximaQ_annuels)
 title('Maxima annuel des débits')
 ylabel('Débit [m3/s]')
 xlabel('Année')
 pause

 hist(maximaQ_annuels,5)
 ylabel('Nombre dannées')
 xlabel('Débit maximal annuel [m3/s]')
 pause

 [parametresQ, parametresQIC] = gevfit(maximaQ_annuels) ;

 parametresQ    % Parametre xi, sigma, mu
 parametresQIC   % intervalles de confiance à 95%
 pause

 dataSortQ = sort(maximaQ_annuels);
 for i = 1:length(dataSortQ)
     FQ(i) = (i-0.5)/(length(dataSortQ));
 end

 stairs(dataSortQ,FQ), hold on
 plot(dataSortQ, gevcdf(dataSortQ,parametresQ(1),parametresQ(2),parametresQ(3)),'g')
 plot(dataSortQ, gevcdf(dataSortQ,parametresQIC(1,1),parametresQIC(1,2),parametresQIC(1,3)),'r')
 plot(dataSortQ, gevcdf(dataSortQ,parametresQIC(2,1),parametresQIC(2,2),parametresQIC(2,3)),'r'), hold off
 title('Fonction de probabilité empirique et estimée')
 xlabel('Débit maximal annuel [m3/s]')
 legend('probabilité empirique','probabilité estimée','intervalle de confiance 95%',4)
 pause

 [parametresGumbel, parametresGumbelIC] = evfit(maximaQ_annuels) ;

 parametresGumbel    % Parametre mu sigma
 parametresGumbelIC   % intervalles de confiance à 95%
 pause

 plot(dataSortQ,dataSortQ), hold on
 plot(gevinv(FQ,parametresQ(1),parametresQ(2),parametresQ(3)),dataSortQ,'.'), hold off
 xlabel('Quantile selon Weibull [m3/s]')
 ylabel('Quantile empirique [m3/s]')
 pause

 tps_rQ = 2:1000;
 tps_r_empiriquesQ = 1./(1-FQ);
 quantile_Q = gevinv((1-1./tps_rQ),parametresQ(1),parametresQ(2),parametresQ(3));
 quantile_Qplus = gevinv((1-1./tps_rQ),parametresQIC(2,1),parametresQIC(2,2),parametresQIC(2,3));
 quantile_Qmoins = gevinv((1-1./tps_rQ),parametresQIC(1,1),parametresQIC(1,2),parametresQIC(1,3));
 
 semilogx(tps_r_empiriquesQ,dataSortQ,'.'), hold on
 semilogx(tps_rQ,quantile_Q)
 plot([100 100],[0 50],'k')
 plot([1 100],[quantile_Q(99) quantile_Q(99)],'k:'), hold off
 axis([1 1000 0 50])
 legend('Quantiles empiriques','Quantiles selon Weibull',2)
 xlabel('Période de retour [an]')
 ylabel('Débits [m3/s]')
 pause

 semilogx(tps_r_empiriquesQ,dataSortQ,'.'), hold on
 semilogx(tps_rQ,quantile_Q)
 semilogx(tps_rQ,quantile_Qplus,'r')
 semilogx(tps_rQ,quantile_Qmoins,'r')
 semilogx(tps_r_empiriquesQ,dataSortQ,'.')
 plot([100 100],[0 60],'k')
 plot([1 100],[quantile_Q(99) quantile_Q(99)],'k:')
 plot([1 100],[quantile_Qplus(99) quantile_Qplus(99)],'k:')
 plot([1 100],[quantile_Qmoins(99) quantile_Qmoins(99)],'k:'), hold off
 axis([1 1000 0 60])
 legend('Quantiles empiriques','Quantiles selon Weibull','Quantiles "limites" selon Weibull 95%',2)
 xlabel('Période de retour [an]')
 ylabel('Débits [m3/s]')
 pause

 [ quantile_Q(99) quantile_Qplus(99) quantile_Qmoins(99)]
 pause

 %%%%%%%%%
 % B - POT
 %%%%%%%%%

 f1 = figure ; hold on ;
 f2 = figure ; hold on ;

 %% Espérance conditionnelle
 for i = 1:100
     s(i)=max(debit)/100 * i;
     Y = debit(find(debit>s(i))) ;
     E_s(i) = mean(Y-s(i)) ;
 end

 % Choix "visuel"
 clear X;
 clear Y;
 X(1) = 0 ; X(2) = 20;
 Y(1) = 7.9 ; Y(2) = 1.9;

 plot(s,E_s,'.k'), hold on
 plot(X,Y,'r'), hold off

 [((Y(2)-Y(1))/(X(2)-X(1))) Y(1)] %xi, sigma

 figure(f1) ; plot(s,E_s,'.k')

 %% Définition du seuil

 choix_seuil = [5 10 15 20]

 for k = 1 : length(choix_seuil)

     seuil = choix_seuil(k)

     clear Y
     Y   = debit(find(debit>seuil)) -seuil ;
     n_s = length(Y) ;
     n_a = length(debit) ;

     [gpparamQ , gpparamQ_ci] = gpfit(Y)

     xi_pareto = gpparamQ(1) ; sigma_pareto = gpparamQ(2) ;

     x(1) = seuil ; x(2) = max(s) ;
     y = (sigma_pareto + xi_pareto * (x-seuil) )/(1-xi_pareto) ;

     figure(f1)
     plot( x , y, '-r') ;
     text( x(2), y(2) , sprintf('seuil %3.1f mm \\rightarrow',seuil),'HorizontalAlignment','right')  ;
      xlabel('seuil') ;
      ylabel('E[X-s|X>s]') ;

     figure(f2)
     lambda = length(debit)/length(Y) / 365;
     plot(     sort(Y+seuil)   ,  1./(1-gpcdf(sort(Y),xi_pareto,sigma_pareto))*lambda       ,'-r') ;
     text( max(sort(Y+seuil))  ,  1./(1-gpcdf(max(sort(Y)),xi_pareto,sigma_pareto))*lambda  ,...
          sprintf('seuil %3.1f mm \\rightarrow',seuil),'HorizontalAlignment','right')
      xlabel('Débit en m3/s') ;
      ylabel('Temps de retour en années') ;
 end
 pause

 % pour la suite nous choisissons un seuil de 30 mm, bien représentatif du
 % "faisseau" 10-40.

     seuil = 15;

     clear Y
     Y   = debit(find(debit>seuil)) -seuil ;
     n_s = length(Y) ;
     n_a = length(debit) ;

     [gpparamQ , gpparamQ_ci] = gpfit(Y)

     xi_pareto = gpparamQ(1) ; sigma_pareto = gpparamQ(2) ;

     x(1) = seuil ; x(2) = max(s) ;
     y = (sigma_pareto + xi_pareto * (x-seuil) )/(1-xi_pareto) ;

     plot(s,E_s,'.k'), hold on
     plot( x , y, '-r'), hold off
      xlabel('seuil') ;
      ylabel('E[X-s|X>s]') ;
     pause

     lambda = length(pluie)/length(Y) / 365;

     plot(sort(Y+seuil),1./(1-gpcdf(sort(Y),xi_pareto,sigma_pareto))*lambda,'-r')
       xlabel('Débits en m3/s') ;
      ylabel('Temps de retour en années') ;
      pause

 dataSort = sort(Y);
 clear F
 for i = 1:length(dataSort)
     F(i) = (i-0.5)/(length(dataSort));
 end

      plot(dataSort+seuil,dataSort+seuil), hold on
 plot(gpinv(F,xi_pareto,sigma_pareto)+seuil,dataSort+seuil ,'.'), hold off
 xlabel('Quantile selon Pareto [m3/s]')
 ylabel('Quantile empirique [m3/s]')
 pause

 % pour la pseudo-période de retour, nous avons :

 m_rQ = 2:10000;
 m_r_empiriquesQ = 1./(1-F);
 quantile_Q_pareto = gpinv((1-1./m_rQ),xi_pareto,sigma_pareto)+seuil;
 pause

 semilogx(m_r_empiriquesQ,dataSort+seuil,'.'), hold on
 semilogx(m_rQ,quantile_Q_pareto), hold off
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Pseudo-période de retour m')
 ylabel('Débits [m3/s]')
 pause

 % On converti cette loi en loi des valeurs extrêmes

 xi = xi_pareto;
 mu = seuil + sigma_pareto/xi_pareto*((n_s/25)^xi_pareto - 1);
 sigma = sigma_pareto * (n_s/25)^xi_pareto;

 xi_plus = gpparamQ_ci(2,1);
 mu_plus = seuil + gpparamQ_ci(2,2)/gpparamQ_ci(2,1)*((n_s/25)^gpparamQ_ci(2,1) - 1);
 sigma_plus = gpparamQ_ci(2,2) * (n_s/25)^gpparamQ_ci(2,1);

 xi_moins = gpparamQ_ci(1,1);
 mu_moins = seuil + gpparamQ_ci(1,2)/gpparamQ_ci(1,1)*((n_s/25)^gpparamQ_ci(1,1) - 1);
 sigma_moins = gpparamQ_ci(1,2) * (n_s/25)^gpparamQ_ci(1,1);

 [xi sigma mu]
 [xi_plus sigma_plus mu_plus]
 [xi_moins sigma_moins mu_moins]
 pause

 tps_rQ = 2:1000;
 quantile_Q_pareto = gevinv((1-1./tps_rQ),xi,sigma,mu);
 quantile_Qplus_pareto = gevinv((1-1./tps_rQ),xi_plus,sigma_plus,mu_plus);
 quantile_Qmoins_pareto = gevinv((1-1./tps_rQ),xi_moins,sigma_moins,mu_moins);

 semilogx(m_r_empiriquesQ*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_rQ,quantile_Q_pareto)
 plot([100 100],[0 50],'k')
 plot([0.1 100],[quantile_Q_pareto(99) quantile_Q_pareto(99)],'k:'), hold off
 axis([0.2 1000 15 50])
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Période de retour [an]')
 ylabel('Débits [m3/s]')
 pause

 semilogx(m_r_empiriquesQ*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_rQ,quantile_Q_pareto)
 semilogx(tps_rQ,quantile_Qplus_pareto,'r')
 semilogx(tps_rQ,quantile_Qmoins_pareto,'r')
 plot([100 100],[0 50],'k')
 plot([0.1 100],[quantile_Q_pareto(99) quantile_Q_pareto(99)],'k:')
 plot([0.1 100],[quantile_Qplus_pareto(99) quantile_Qplus_pareto(99)],'k:')
 plot([0.1 100],[quantile_Qmoins_pareto(99) quantile_Qmoins_pareto(99)],'k:'), hold off
 axis([0.2 1000 15 50])
 legend('Quantiles empiriques','Quantiles selon Pareto',2)
 xlabel('Période de retour [an]')
 ylabel('Débits [m3/s]')
 pause

 [ quantile_Q_pareto(99) quantile_Qplus_pareto(99) quantile_Qmoins_pareto(99)]
 pause

 %Comparaison Frechet - Pareto :

 semilogx(m_r_empiriquesQ*lambda,dataSort+seuil,'.'), hold on
 semilogx(tps_rQ,quantile_Q_pareto)
 semilogx(tps_rQ,quantile_Q,'r')
 plot([100 100],[15 50],'k')
 plot([0.1 100],[quantile_Q(99) quantile_Q(99)],'k:')
 plot([0.1 100],[quantile_Q_pareto(99) quantile_Q_pareto(99)],'k:'), hold off
 axis([0.2 1000 20 40])
 legend('Quantiles empiriques','Quantiles selon Pareto','Quantiles selon Weibull',2)
 xlabel('Période de retour [an]')
 ylabel('Débit [m3/s]')

 %%%%%%%%%
 % C
 %%%%%%%%%

 % bla bla

 %%%%%%%%%
 % D
 %%%%%%%%%

 plot(debit,pluie,'.k')
 xlabel('Débits [m3/s]')
 ylabel('Pluies [mm]')
 title('Couples pluie/débit pour la chronique de 1974 à 1998')