source: trunk/STATISTICS/c_timecorrelate.pro @ 2

Last change on this file since 2 was 2, checked in by opalod, 22 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.0 KB
Line 
1;+
2; NAME:
3;       C_TIMECORRELATE
4;
5; PURPOSE:
6;       This function computes the "time cross correlation" Pxy(L) or
7;       the "time cross covariance" between 2 arrays (this is some
8;       kind of c_correlate but for multidimenstionals arrays) as a
9;       function of the lag (L).
10;
11; CATEGORY:
12;       Statistics.
13;
14; CALLING SEQUENCE:
15;       Result = c_timecorrelate(X, Y, Lag)
16;
17; INPUTS:
18;       X:   an Array which last dimension is the time dimension of
19;       size n, float or double.
20;
21;       Y:    an Array which last dimension is the time dimension of
22;       size n, float or double.
23;
24;     LAG:    A scalar or n-element vector, in the interval [-(n-2), (n-2)],
25;             of type integer that specifies the absolute distance(s) between
26;             indexed elements of X.
27;
28; KEYWORD PARAMETERS:
29;       COVARIANCE:    If set to a non-zero value, the sample cross
30;                      covariance is computed.
31;
32;       DOUBLE:        If set to a non-zero value, computations are done in
33;                      double precision arithmetic.
34;
35; EXAMPLE
36;       Define two n-element sample populations.
37;         x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
38;         y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95]
39;
40;       Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7
41;         lag = [-5, 0, 1, 5, 6, 7]
42;         result = c_timecorrelate(x, y, lag)
43;
44;       The result should be:
45;         [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
46;
47; PROCEDURE:
48;
49;
50;   FOR L>=0
51;
52;                              n-L-1
53;                              sigma  (X[k]-Xmean)(Y[k+L]-Ymean)
54;                               k=0
55;   correlation(X,Y,L)=------------------------------------------------------
56;                              n-1                     n-1
57;                      sqrt( (sigma  (X[k]-Xmean)^2)*(sigma  (Y[k]-Ymean)^2))
58;                              k=0                     k=0
59;
60;
61;
62;                              n-L-1
63;                              sigma  (X[k]-Xmean)(Y[k+L]-Ymean)
64;                               k=0
65;   covariance(X,Y,L)=------------------------------------------------------
66;                                            n
67;
68;   FOR L<0
69;
70;
71;                              n-L-1
72;                              sigma  (X[k+L]-Xmean)(Y[k]-Ymean)
73;                               k=0
74;   correlation(X,Y,L)=------------------------------------------------------
75;                              n-1                     n-1
76;                      sqrt( (sigma  (X[k]-Xmean)^2)*(sigma  (Y[k]-Ymean)^2))
77;                              k=0                     k=0
78;
79;
80;
81;                              n-L-1
82;                              sigma  (X[k+L]-Xmean)(Y[k]-Ymean)
83;                               k=0
84;   covariance(X,Y,L)=------------------------------------------------------
85;                                            n
86;
87;   Where Xmean and Ymean are the time means of the sample populations
88;   x=(x[t=0],x[t=1],...,x[t=n-1]) and y=(y[t=0],y[t=1],...,y[t=n-1]),
89;   respectively.
90;
91;
92;
93; REFERENCE:
94;       INTRODUCTION TO STATISTICAL TIME SERIES
95;       Wayne A. Fuller
96;       ISBN 0-471-28715-6
97;
98; MODIFICATION HISTORY:
99;       01/03/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
100;       Based on the C_CORRELATE procedure of IDL
101;-
102
103FUNCTION TimeCross_Cov, X, Y, M, nT, Double = Double, ZERO2NAN = zero2nan
104;
105   if double then one = 1.0d ELSE one = 1.0
106; Sample cross covariance function.
107   TimeDim = size(X, /n_dimensions)
108   Xmean = TOTAL(X, TimeDim, Double = Double) / nT
109   Xmean = Xmean[*]#replicate(one,  nT - M)
110   Ymean = TOTAL(Y, TimeDim, Double = Double) / nT
111   Ymean = Ymean[*]#replicate(one,  nT - M)
112;
113   case TimeDim of
114      1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (Y[M:nT - 1L] - Ymean) $
115                    , Double = Double)
116      2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * (Y[*, M:nT - 1L] - Ymean) $
117                    , TimeDim, Double = Double)
118      3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, M:nT - 1L] - Ymean) $
119                    , TimeDim, Double = Double)
120      4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, *, M:nT - 1L] - Ymean) $
121                    , TimeDim, Double = Double)
122   ENDCASE
123   if keyword_set(zero2nan) then begin
124      zero = where(res EQ 0)
125      if zero[0] NE -1 then res[zero] = !values.f_nan
126   ENDIF
127;
128   RETURN, res
129
130END
131
132FUNCTION C_Timecorrelate, X, Y, Lag, Covariance = Covariance, Double = Double
133
134;Compute the sample cross correlation or cross covariance of
135;(Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l).
136
137   ON_ERROR, 2
138
139   xsize = SIZE(X)
140   ysize = SIZE(Y)
141   nt = float(xsize[xsize[0]])
142   NDim = xsize[0]
143
144   if total(xsize[0:xsize[0]] NE ysize[0:ysize[0]]) NE 0 then $
145    MESSAGE, "X and Y arrays must have the same size and the same dimensions"
146
147;Check length.
148   if nt lt 2 then $
149    MESSAGE, "Time dimension of X and Y arrays must contain 2 or more elements."
150   
151;If the DOUBLE keyword is not set then the internal precision and
152;result are identical to the type of input.
153   if N_ELEMENTS(Double) eq 0 then $
154    Double = (Xsize[Xsize[0]+1] eq 5 or ysize[ysize[0]+1] eq 5)
155
156   if n_elements(lag) EQ 0 then lag = 0
157   nLag = N_ELEMENTS(Lag)
158
159   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
160
161   case NDim of
162      1:if Double eq 0 then Correl = FLTARR(nLag) else Correl = DBLARR(nLag)
163      2:if Double eq 0 then Correl = FLTARR(Xsize[1], nLag) else Correl = DBLARR(Xsize[1], nLag)
164      3:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], nLag) $
165      else Correl = DBLARR(Xsize[1], Xsize[2], nLag)
166      4:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $
167      else Correl = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag)
168   endcase
169
170   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross Correlation.
171      for k = 0, nLag-1 do begin
172         if Lag[k] ge 0 then BEGIN
173            case NDim of
174               1:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
175                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
176                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
177               2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
178                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
179                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
180               3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
181                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
182                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
183               4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
184                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
185                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
186            endcase
187         ENDIF else BEGIN
188            case NDim of
189               1:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
190                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
191                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
192               2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
193                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
194                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
195               3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
196                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
197                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
198               4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
199                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
200                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
201            endcase
202         ENDELSE
203      endfor
204   endif else begin             ;Compute Cross Covariance.
205      for k = 0, nLag-1 do begin
206         if Lag[k] ge 0 then BEGIN
207            case NDim of
208               1:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
209               2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
210               3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
211               4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
212            ENDCASE
213         ENDIF else BEGIN
214            case NDim of
215               1:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
216               2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
217               3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
218               4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
219            ENDCASE
220         ENDELSE
221      endfor
222   endelse
223
224   if Double eq 0 then RETURN, FLOAT(Correl) else $
225    RETURN, Correl
226
227END
228 
Note: See TracBrowser for help on using the repository browser.