;+ ; 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