source: trunk/RESULTATS/regandsignif_seriestemporelles.m @ 29

%REGANDSIGNIF_SERIESTEMPORELLES calcule la regression lineaire d'une serie temporelle issue du champ de donnees initial
% (selectionne en un point d'espace particulier) et une composante principale choisie      
%
% Date: Avril-Mai 2008
% Auteur: juliette.mignot@locean-ipsl.upmc.fr
% Contexte Ecole d'ete ENVI-STAT 2008
%+
%
% module
% ======
%
% ``regandsignif_seriestemporelles``
%
% DESCRIPTION
% ===========
%
% ``regandsignif_seriestemporelles`` lit les fichiers fichiers 
% ``../DONNEES/lsmask_tropatl_PacandMedblancs.mat``,
% ``../DONNEES/skt.mon.tropatl.mat``.
%
% ``regandsignif_seriestemporelles`` lit aussi le fichier 
% ``../DONNEES//eof_SSTan_tropAtl30N-20S.mat``
% produit par eof_NCEP_an.m_.
%
% ``regandsignif_seriestemporelles`` calcule ++
%
% ``regandsignif_seriestemporelles`` affiche ++
%
% ``regandsignif_seriestemporelles`` sauve cette image dans 
% ``./eof_SSTan_tropAtl30N-20S.mat_regandsignif_seriestemporelles_1.ps``.
%
% EXAMPLES
% ========
%
% ::
%
% >> addpath('../PROGRAMMES/')
% >> tpacpandreg_startup
% >> regandsignif_seriestemporelles
%
%
% SEE ALSO
% ========
%
% tpacpandreg_startup.m_
%
% .. _tpacpandreg_startup.m : tpacpandreg_startup.m.html
%
% eof_NCEP_an.m_
%
% .. _eof_NCEP_an.m : eof_NCEP_an.m.html
%
% TODO
% ====
%
% improve description
%
% pourquoi ``return`` à la fin
%
% ajout xlabel et ylabel
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% 
% - fplod 2009-08-26T11:00:35Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * remplacement de ``eps`` par ``ps`` pour pouvoir imprimer les figures 
%     produites par ``octave``
% 
% - fplod 2009-08-25T10:40:16Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * modif commentaires liée au renommage
%
% - jmignot 2009-08-25
%
%   * RESULTATS/Reg_seriestemporelles.m renommé en 
%     RESULTATS/regandsignif_seriestemporelles.m pour éviter confusion
%     avec PROGRAMMES/Reg_seriestemporelles.m
%   * chgt de noms de fichiers : ``EOF_`` devient ``eof_``
%
% - fplod 2009-08-24T11:44:38Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add header
%   * répertoire de données relatifs au répertoire courant
%   * suppression de ``clear all`` présent dans ``tpacpandreg_startup.m``
%   * as octave save default format is ASCII and no default extension
%     while matlab save default format is MAT v5 mat-file (little endian) and
%     default extension is ``.mat``, write more precise save instruction
%   * save the figure
%   * replace file radical ``EOF_SST_tropAtl`` by ``EOF_SSTan_tropAtl``
%   * split lines with multiple statements
%   * handle ``corrcoef`` octave/matlab difference of output
%
%-

%%
%% PARAMETRES DU PROGRAMME
% selection de la composante principale a considerer
ipc=1; %rang de la composante principale que l'on va considerer ici

% les coordonnees du point d'espace auquel on veut comparer cette
% composante principale doit etre choisi plus bas, apres chargement des
% longitudes et latitudes
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%              CHARGEMENT ET PREPARATION DES DONNEES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nannees=60;
repertoire=[ '..' filesep 'DONNEES' filesep ]; % a adapter a chaque cas: chemin d acces aux donnees

% Mask
%------
file='lsmask_tropatl_PacandMedblancs.mat'
   load([repertoire file]);
% Data
%--------
  repertoire2=[repertoire 'MAT' filesep] % a adapter a chaque cas: chemin d acces aux donnees
  file=['skt.mon.tropatl.mat']
  load([repertoire2 file]); % chargement des donnees
  tab=permute(tab,[2 3 1]); % permutation des dimensions pour organiser le tableau sous la forme [lat lon temps]

% Application du masque
%----------------------
mask3D=repmat(mask,[1 1 size(tab,3)]);
tab(mask3D~=0)=NaN;


%  MANIPULATIONS TEMPORELLES: MOYENNES 

% Construction de moyennes annuelles
%---------------------------------------
tab=reshape(tab,size(tab,1),size(tab,2),12,nannees);
tab_an=squeeze(mean(tab,3));
%Calcul des anomalies
%tab_an_mean=squeeze(mean(tab_an,3));
%tab_an_mean=repmat(tab_an_mean,[1 1 size(tab_an,3)]);
%tab_an=tab_an-tab_an_mean;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% selection d'un point d'espace
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% le point d'espace est selectionne en reperant son indice i,j dans le
% tableau de longitude et de latitude

%i=find(lon>-96 & lon<-94);
%j=find(lat>19 & lat<21);
i=find(lon>-96 & lon<-94);
j=find(lat>19 & lat<21);
%i=find(lon>-21 & lon<-19);
%j=find(lat>19 & lat<21);
%i=find(lon>4 & lon<6)
%j=find(lat>-11 & lat<-9)

SST=squeeze(tab_an(j,i,:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% chargement des composantes principales
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

file=['eof_SSTan_tropAtl30N-20S.mat'];
titledataeof='SST' % je definis ce nom pour pouvoir l'inclure dans le graphique plus bas. 
if (run_octave == 0)
 load([repertoire file]);
else
 load('-v7',[repertoire file]);
end
PC=PC(:,ipc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%               NUAGE DE POINTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ifigure=1;
figure(ifigure); 
plot(SST,PC,'x')
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%   REGRESSION LINEAIRE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%[b0,b1,s,R,sigb0,sigb1] = linreg(SST,serie);

SSTmean=mean(SST);PCmean=mean(PC);
n=length(SST);

% calcul de b1 et b0, coefficient directeur et ordonnee a l origine de la
% droite de regression
b1=sum((PC-PCmean).*(SST-SSTmean))/(sum(SST.*SST)-n*SSTmean^2);
b0=PCmean-b1*SSTmean;

PCc=b0+b1*SST;

%trace de la droite de regression
%---------------------------------
plot(SST,PCc,'r')

%Calcul du coefficient de détermination
%----------------------------------------
s=sqrt(sum((PC-PCc).^2)/(n-2));
R=sum((PCc-PCmean).^2)/sum((PC-PCmean).^2);

%significativite statistique comparaison de la correlation r et de la
%valeur critique rc donnee par la table de student en fonction du nombre de
%degres de liberte ndf
%------------------------------

r=corrcoef(SST,PC);
if (run_octave == 0)
 r=r(2,1);
end

% remarque: b1=r*std(PC)/std(SST)

% on estime le nombre de degrés de liberté
ndf=nannees; %(l'atmosphère est approximativement un bruit blanc à l'échelle de temps interanuellee)

% on calcule le seuil de significativité au percentile P en inversant
% l'écart réduit T=y*sqrt(nn)./(1-y.^2).^0.5 qui est distribue selon
% un loi de student avec nn=ndf-2 degres de liberte.
P=.95; 
nn=ndf-2;
t=tinv(P,nn);
rc = t ./ ( nn + t.^2 ).^0.5;
% rc est la valeur ccritique de la correlation pour rejeter l'hypothese nulle

title(['r = ' num2str(r) '  valeur critique: rc= ' num2str(rc)])

% significativite statistique de b0 et b1
%--------------------------------------------

% sauvegarde de la figure
printer='ps';
print_printer=['-d', printer];
fullfilename=['.' filesep file '_' mfilename '_'  num2str(ifigure) '.' printer];
print(print_printer,fullfilename);
clear print_printer;
clear fullfilename;
clear printer;
%
clear ifigure;

return