correl.pro (Kronos IDL library)

;+
; Contains the correl procedure
;
; :Author:
;    Laurent Lamy or Philippe Zarka
;
; :History:
;    2009/11/24: Created
;
;    2009/11/24: Last Edit
;-

;+
; Calcule le coefficient de correlation C et son degre 
; de confiance P de 2 tableaux X(N) et Y(N).
;       $$C = \frac{\sigma_{XY}}{\sigma_{X} * \sigma_{Y}}$$
; avec $$\sigma_X^{2} = \frac{1}{N}  \sum_{i=1}^{N} (X_{i}-X_{m})^{2}$$
;      $$\sigma_{Y}^{2} = \frac{1}{N}  \sum_{i=1}^{N} (Y_{i}-Y_{m})^{2}$$
;      $$\sigma_{XY} = \frac{1}{N}  \sum_{i=1}^{N} (X_{i}-X_{m})*(Y_{i}-Y_{m})$$
; et   $$P = 100 * \left|\frac{2}{\sqrt{\pi}}  \frac{\Gamma(\frac{N-1}{2})}{\Gamma(\frac{N-2}{2})} \int_0^C (1-u^{2})^{\frac{N-4}{2}} \mathrm{d}u \right|$$
;
; Pour le calcul de P, on opere de la maniere suivante:
;
; 1) $$\frac{\Gamma(\frac{N-1}{2})}{\Gamma(\frac{N-2}{2})}$$  est calculé avec la fonction 
;   LNGAMMA (erreur < 1%)
;
; 2)  $$\int_0^C (1-u^{2})^{\frac{N-4}{2}} \mathrm{d}u $$ est calculée par recurrence:
;   si $$I_{k} = \int_0^C (1-u^{2})^{k} \mathrm{d}u $$
;   $$ I_{k} = \frac{1}{2k+1} C(1-C^{2})^{k} + \frac{2k}{2k+1} I_{k-1} $$
;   $$I_{0} = C$$
;   $$I_{1/2} = \frac{1}{2} \arcsin C + \frac{C}{2} \sqrt{1-C^{2}}$$
;
; Le calcul de P n'est pas effectue pour N > 23000.
;
; :Params:
;    X: in, required, type=fltarr
;       Tableau
;    Y: in, required, type=fltarr
;       Tableau
;    C: out, required, type=float
;       coefficient de correlation (-1.<=C<=1)
;    P: out, required, type=float
;       degre de confiance ( 0.<=P<=100.)
;-
pro CORREL, X, Y, C, P
  on_error,1
  C=0.	& P=999.
  N=n_elements(X)
  Xm=float(total(X)/N)		& Ym=float(total(Y)/N)
  sigma_X=sqrt(total((X-Xm)^2)/N)	& sigma_Y=sqrt(total((Y-Ym)^2)/N)
  sigma_XY=total((X-Xm)*(Y-Ym))/N
  if sigma_X eq 0. or sigma_Y eq 0. then return
  C=float(sigma_XY/sigma_X/sigma_Y)
;	print,n,xm,ym,sigma_x,sigma_y,sigma_xy,c
  if N gt 100000 then return
  R=2/sqrt(!Pi)*exp(lngamma(double(N-1)/2.)-lngamma(double(N-2)/2.))
  if (N mod 2) eq 0 then begin
	ii=double(C)
	k=1.
  endif
  if (N mod 2) eq 1 then begin
	ii=.5*asin(double(C))+C/2.d0*sqrt(1.d0-C*C)
	k=1.5
  endif
  i0=1000.
  while (k le (N-4.)/2.) and (abs(ii-i0)/i0 ge 1.d-5) do begin
	   i0=ii
	   ii=(C*(1-C*C)^k+2.*k*i0)/(2.*k+1.)
	   k=k+1.
  endwhile
  P=100.*abs(R*ii)
;	print,'c,p,k,ii=',c,p,k,ii
return
end