1 | %REGANDSIGNIF_SERIESTEMPORELLES calcule la regression lineaire d'une serie temporelle issue du champ de donnees initial |
---|
2 | % (selectionne en un point d'espace particulier) et une composante principale choisie |
---|
3 | % |
---|
4 | % Date: Avril-Mai 2008 |
---|
5 | % Auteur: juliette.mignot@locean-ipsl.upmc.fr |
---|
6 | % Contexte Ecole d'ete ENVI-STAT 2008 |
---|
7 | %+ |
---|
8 | % |
---|
9 | % module |
---|
10 | % ====== |
---|
11 | % |
---|
12 | % ``regandsignif_seriestemporelles`` |
---|
13 | % |
---|
14 | % DESCRIPTION |
---|
15 | % =========== |
---|
16 | % |
---|
17 | % ``regandsignif_seriestemporelles`` lit les fichiers fichiers |
---|
18 | % ``../DONNEES/lsmask_tropatl_PacandMedblancs.mat``, |
---|
19 | % ``../DONNEES/skt.mon.tropatl.mat``. |
---|
20 | % |
---|
21 | % ``regandsignif_seriestemporelles`` lit aussi le fichier |
---|
22 | % ``../DONNEES//eof_SSTan_tropAtl30N-20S.mat`` |
---|
23 | % produit par eof_NCEP_an.m_. |
---|
24 | % |
---|
25 | % ``regandsignif_seriestemporelles`` calcule ++ |
---|
26 | % |
---|
27 | % ``regandsignif_seriestemporelles`` affiche ++ |
---|
28 | % |
---|
29 | % ``regandsignif_seriestemporelles`` sauve cette image dans |
---|
30 | % ``./eof_SSTan_tropAtl30N-20S.mat_regandsignif_seriestemporelles_1.ps``. |
---|
31 | % |
---|
32 | % EXAMPLES |
---|
33 | % ======== |
---|
34 | % |
---|
35 | % :: |
---|
36 | % |
---|
37 | % >> addpath('../PROGRAMMES/') |
---|
38 | % >> tpacpandreg_startup |
---|
39 | % >> regandsignif_seriestemporelles |
---|
40 | % |
---|
41 | % |
---|
42 | % SEE ALSO |
---|
43 | % ======== |
---|
44 | % |
---|
45 | % tpacpandreg_startup.m_ |
---|
46 | % |
---|
47 | % .. _tpacpandreg_startup.m : tpacpandreg_startup.m.html |
---|
48 | % |
---|
49 | % eof_NCEP_an.m_ |
---|
50 | % |
---|
51 | % .. _eof_NCEP_an.m : eof_NCEP_an.m.html |
---|
52 | % |
---|
53 | % TODO |
---|
54 | % ==== |
---|
55 | % |
---|
56 | % improve description |
---|
57 | % |
---|
58 | % pourquoi ``return`` à la fin |
---|
59 | % |
---|
60 | % ajout xlabel et ylabel |
---|
61 | % |
---|
62 | % EVOLUTIONS |
---|
63 | % ========== |
---|
64 | % |
---|
65 | % $Id$ |
---|
66 | % |
---|
67 | % - fplod 2009-09-14T14:10:07Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
68 | % |
---|
69 | % * explicit MarkerSize value to avoid difference between |
---|
70 | % octave 3.0.2 and octave 3.2.2 plots |
---|
71 | % * explicit X and Y limits to avoid difference between |
---|
72 | % matlab, octave 3.0.2 and octave 3.2.2 plots |
---|
73 | % |
---|
74 | % - fplod 2009-08-26T11:00:35Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
75 | % |
---|
76 | % * remplacement de ``eps`` par ``ps`` pour pouvoir imprimer les figures |
---|
77 | % produites par ``octave`` |
---|
78 | % |
---|
79 | % - fplod 2009-08-25T10:40:16Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
80 | % |
---|
81 | % * modif commentaires liée au renommage |
---|
82 | % |
---|
83 | % - jmignot 2009-08-25 |
---|
84 | % |
---|
85 | % * RESULTATS/Reg_seriestemporelles.m renommé en |
---|
86 | % RESULTATS/regandsignif_seriestemporelles.m pour éviter confusion |
---|
87 | % avec PROGRAMMES/Reg_seriestemporelles.m |
---|
88 | % * chgt de noms de fichiers : ``EOF_`` devient ``eof_`` |
---|
89 | % |
---|
90 | % - fplod 2009-08-24T11:44:38Z aedon.locean-ipsl.upmc.fr (Darwin) |
---|
91 | % |
---|
92 | % * add header |
---|
93 | % * répertoire de données relatifs au répertoire courant |
---|
94 | % * suppression de ``clear all`` présent dans ``tpacpandreg_startup.m`` |
---|
95 | % * as octave save default format is ASCII and no default extension |
---|
96 | % while matlab save default format is MAT v5 mat-file (little endian) and |
---|
97 | % default extension is ``.mat``, write more precise save instruction |
---|
98 | % * save the figure |
---|
99 | % * replace file radical ``EOF_SST_tropAtl`` by ``EOF_SSTan_tropAtl`` |
---|
100 | % * split lines with multiple statements |
---|
101 | % * handle ``corrcoef`` octave/matlab difference of output |
---|
102 | % |
---|
103 | %- |
---|
104 | |
---|
105 | %% |
---|
106 | %% PARAMETRES DU PROGRAMME |
---|
107 | % selection de la composante principale a considerer |
---|
108 | ipc=1; %rang de la composante principale que l'on va considerer ici |
---|
109 | |
---|
110 | % les coordonnees du point d'espace auquel on veut comparer cette |
---|
111 | % composante principale doit etre choisi plus bas, apres chargement des |
---|
112 | % longitudes et latitudes |
---|
113 | %% |
---|
114 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
115 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
116 | %%%% CHARGEMENT ET PREPARATION DES DONNEES |
---|
117 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
118 | |
---|
119 | nannees=60; |
---|
120 | repertoire=[ '..' filesep 'DONNEES' filesep ]; % a adapter a chaque cas: chemin d acces aux donnees |
---|
121 | |
---|
122 | % Mask |
---|
123 | %------ |
---|
124 | file='lsmask_tropatl_PacandMedblancs.mat' |
---|
125 | load([repertoire file]); |
---|
126 | % Data |
---|
127 | %-------- |
---|
128 | repertoire2=[repertoire 'MAT' filesep] % a adapter a chaque cas: chemin d acces aux donnees |
---|
129 | file=['skt.mon.tropatl.mat'] |
---|
130 | load([repertoire2 file]); % chargement des donnees |
---|
131 | tab=permute(tab,[2 3 1]); % permutation des dimensions pour organiser le tableau sous la forme [lat lon temps] |
---|
132 | |
---|
133 | % Application du masque |
---|
134 | %---------------------- |
---|
135 | mask3D=repmat(mask,[1 1 size(tab,3)]); |
---|
136 | tab(mask3D~=0)=NaN; |
---|
137 | |
---|
138 | |
---|
139 | % MANIPULATIONS TEMPORELLES: MOYENNES |
---|
140 | |
---|
141 | % Construction de moyennes annuelles |
---|
142 | %--------------------------------------- |
---|
143 | tab=reshape(tab,size(tab,1),size(tab,2),12,nannees); |
---|
144 | tab_an=squeeze(mean(tab,3)); |
---|
145 | %Calcul des anomalies |
---|
146 | %tab_an_mean=squeeze(mean(tab_an,3)); |
---|
147 | %tab_an_mean=repmat(tab_an_mean,[1 1 size(tab_an,3)]); |
---|
148 | %tab_an=tab_an-tab_an_mean; |
---|
149 | |
---|
150 | %% |
---|
151 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
152 | %%%% selection d'un point d'espace |
---|
153 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
154 | % le point d'espace est selectionne en reperant son indice i,j dans le |
---|
155 | % tableau de longitude et de latitude |
---|
156 | |
---|
157 | %i=find(lon>-96 & lon<-94); |
---|
158 | %j=find(lat>19 & lat<21); |
---|
159 | i=find(lon>-96 & lon<-94); |
---|
160 | j=find(lat>19 & lat<21); |
---|
161 | %i=find(lon>-21 & lon<-19); |
---|
162 | %j=find(lat>19 & lat<21); |
---|
163 | %i=find(lon>4 & lon<6) |
---|
164 | %j=find(lat>-11 & lat<-9) |
---|
165 | |
---|
166 | SST=squeeze(tab_an(j,i,:)); |
---|
167 | |
---|
168 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
169 | %%%%% chargement des composantes principales |
---|
170 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
171 | |
---|
172 | file=['eof_SSTan_tropAtl30N-20S.mat']; |
---|
173 | titledataeof='SST' % je definis ce nom pour pouvoir l'inclure dans le graphique plus bas. |
---|
174 | if (run_octave == 0) |
---|
175 | load([repertoire file]); |
---|
176 | else |
---|
177 | load('-v7',[repertoire file]); |
---|
178 | end |
---|
179 | PC=PC(:,ipc); |
---|
180 | |
---|
181 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
182 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
183 | %%%% NUAGE DE POINTS |
---|
184 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
185 | |
---|
186 | ifigure=1; |
---|
187 | figure(ifigure); |
---|
188 | plot(SST,PC,'x','markersize', 6) |
---|
189 | xlim([24 26]); |
---|
190 | ylim([-3 3]); |
---|
191 | hold on |
---|
192 | |
---|
193 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
194 | %%%%%% REGRESSION LINEAIRE |
---|
195 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
196 | |
---|
197 | %[b0,b1,s,R,sigb0,sigb1] = linreg(SST,serie); |
---|
198 | |
---|
199 | SSTmean=mean(SST);PCmean=mean(PC); |
---|
200 | n=length(SST); |
---|
201 | |
---|
202 | % calcul de b1 et b0, coefficient directeur et ordonnee a l origine de la |
---|
203 | % droite de regression |
---|
204 | b1=sum((PC-PCmean).*(SST-SSTmean))/(sum(SST.*SST)-n*SSTmean^2); |
---|
205 | b0=PCmean-b1*SSTmean; |
---|
206 | |
---|
207 | PCc=b0+b1*SST; |
---|
208 | |
---|
209 | %trace de la droite de regression |
---|
210 | %--------------------------------- |
---|
211 | plot(SST,PCc,'r') |
---|
212 | |
---|
213 | %Calcul du coefficient de détermination |
---|
214 | %---------------------------------------- |
---|
215 | s=sqrt(sum((PC-PCc).^2)/(n-2)); |
---|
216 | R=sum((PCc-PCmean).^2)/sum((PC-PCmean).^2); |
---|
217 | |
---|
218 | %significativite statistique comparaison de la correlation r et de la |
---|
219 | %valeur critique rc donnee par la table de student en fonction du nombre de |
---|
220 | %degres de liberte ndf |
---|
221 | %------------------------------ |
---|
222 | |
---|
223 | r=corrcoef(SST,PC); |
---|
224 | if (run_octave == 0) |
---|
225 | r=r(2,1); |
---|
226 | end |
---|
227 | |
---|
228 | % remarque: b1=r*std(PC)/std(SST) |
---|
229 | |
---|
230 | % on estime le nombre de degrés de liberté |
---|
231 | ndf=nannees; %(l'atmosphère est approximativement un bruit blanc à l'échelle de temps interanuellee) |
---|
232 | |
---|
233 | % on calcule le seuil de significativité au percentile P en inversant |
---|
234 | % l'écart réduit T=y*sqrt(nn)./(1-y.^2).^0.5 qui est distribue selon |
---|
235 | % un loi de student avec nn=ndf-2 degres de liberte. |
---|
236 | P=.95; |
---|
237 | nn=ndf-2; |
---|
238 | t=tinv(P,nn); |
---|
239 | rc = t ./ ( nn + t.^2 ).^0.5; |
---|
240 | % rc est la valeur ccritique de la correlation pour rejeter l'hypothese nulle |
---|
241 | |
---|
242 | title(['r = ' num2str(r) ' valeur critique: rc= ' num2str(rc)]) |
---|
243 | |
---|
244 | % significativite statistique de b0 et b1 |
---|
245 | %-------------------------------------------- |
---|
246 | |
---|
247 | % sauvegarde de la figure |
---|
248 | printer='ps'; |
---|
249 | print_printer=['-d', printer]; |
---|
250 | fullfilename=['.' filesep file '_' mfilename '_' num2str(ifigure) '.' printer]; |
---|
251 | print(print_printer,fullfilename); |
---|
252 | clear print_printer; |
---|
253 | clear fullfilename; |
---|
254 | clear printer; |
---|
255 | % |
---|
256 | clear ifigure; |
---|
257 | |
---|
258 | return |
---|
259 | |
---|