about CCA

From: Stergios Misios <stergios.misios_at_nyahnyahspammersnyahnyah>
Date: Thu Oct 22 2009 - 12:16:06 MDT

hi all,
i try to perform canonical correlation analysis with the cancor routine.
what i am doing is to read the two fields, to standardize to reshape
from (time,lev,lat) to (time,lev*lat) and finally to apply cancor.
briefly a code snippet,
==============================================
 ;standarize
;----------
print("standarize matrices...")
pw_std=dim_standardize_n(pw,1,0)
uw_std=dim_standardize_n(uwind,1,0)

;reshape arrays
;--------------
print("Reshape matrices")
do itime=0,ntime_pw-1
  trans_pw(:,itime)=ndtooned(pw_std(itime,:,:))
  trans_uw(:,itime)=ndtooned(uw_std(itime,:,:))
end do

;release memory
;-------------
delete(pw_std)
delete(uw_std)

; compute CCA vector
;-------------------
print("Predictor XBs:"+nlev_pw*nlat_pw)
print("Predictand YBs:"+nlev_uw*nlat_uw)
print("apply CCA method...")

opt=False
ccam=cancor(trans_pw,trans_uw,opt)
print("End...")
print(ccam)
=========================================
the output i must say is very strange,

Variable: ccam
Type: float
Total Size: 816 bytes
            204 values
Number of Dimensions: 1
Dimensions and sizes: [204]
Coordinates:
Number Of Attributes: 5
  ndof : <ARRAY of 204 elements>
  coefy : <ARRAY of 41616 elements>
  coefx : <ARRAY of 48960 elements>
  wlam : <ARRAY of 204 elements>
  chisq : <ARRAY of 204 elements>
(0) 111823
(1) 87949.13
(2) 62144.38
(3) 50775.16
(4) 40324.11
(5) 36114.61
(6) 27940.67
(7) 25172.04
(8) 22960.51
(9) 19126.69
(10) 18376.85
(11) 16603.9
(12) 15086.54
(13) 14262.58
(14) 12602.33

This i guess it is suppoded to be the correlation coefficients and must
be <1. where am i wrong?
i guess also i must weight the values but i left it for later...
please find attached the script...

Thank you in advance,
stergios

;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
;------------------------------------------------
exp_name="c2021_blz"
file_pwn_mm="c2021_blz_w)2140-2169_nac_DJF_smm.nc"
file_uw_mm="c2021_blz_u_2140-2169_nac_zm_DJF_smm.nc"
var_wn="wn1"
;------------------------------------------------

;read first file
;------------

pwn_ff=addfile(file_pwn_mm,"r")
pw=pwn_ff->wn1(:,{50:1},{60:80},0)

dsizes_pw=dimsizes(pw)

ntime_pw=dsizes_pw(0)
nlev_pw=dsizes_pw(1)
nlat_pw=dsizes_pw(2)

;read second
;----------
print("read zonal wind")
uw_ff=addfile(file_uw_mm,"r")
uwind=uw_ff->u(:,{3000:100},{-10:10},0)

dsizes_uw=dimsizes(uwind)
;print(dsizes_uw)
ntime_uw=dsizes_uw(0)
nlev_uw=dsizes_uw(1)
nlat_uw=dsizes_uw(2)
;printVarSummary(pw)

;DEFINITIONS
;-----------
;transformed timexlatxlev grid to timexlev*lat
;-------------------------------------------
trans_pw=new((/nlev_pw*nlat_pw,ntime_pw/),float)
trans_uw=new((/nlev_uw*nlat_uw,ntime_uw/),float)

;standarize
;----------
print("standarize matrices...")
pw_std=dim_standardize_n(pw,1,0)
uw_std=dim_standardize_n(uwind,1,0)

;reshape arrays
;--------------
print("Reshape matrices")
do itime=0,ntime_pw-1
  trans_pw(:,itime)=ndtooned(pw_std(itime,:,:))
  trans_uw(:,itime)=ndtooned(uw_std(itime,:,:))
end do

;release memory
;-------------
delete(pw_std)
delete(uw_std)

; compute CCA vector
;-------------------
print("Predictor XBs:"+nlev_pw*nlat_pw)
print("Predictand YBs:"+nlev_uw*nlat_uw)
print("apply CCA method...")

opt=False
ccam=cancor(trans_pw,trans_uw,opt)
print("End...")
print(ccam)
print(ccam@ndof)
print(nlev_uw*nlat_uw)

;reshaped vectors to levxlat grid for plotting
;-------------------------------
cca_pw=new((/nlev_uw*nlat_uw,nlev_pw*nlat_pw/),float)
cca_uw=new((/nlev_uw*nlat_uw,nlev_uw*nlat_uw/),float)
rcca_pw=new((/nlev_pw,nlat_pw/),float)
rcca_uw=new((/nlev_uw,nlat_uw/),float)

;Propability
;-----------
;print("Calculate propability")
;prob = gammainc( 0.5*ccam@chisq, 0.5*ccam@ndof )

cca_pw=ccam@coefx
cca_uw=ccam@coefy
;print(ccam@coefx)

;get levels
;----------
lev_uw=(/uwind&lev/)
lev_pw=(/pw&lev/)

;prepare plot
;------------
plot=new(2,graphic)

;first figure
;------------

res_pw = True

res_pw@gsnDraw = False
res_pw@gsnFrame = False
res_pw@sfXArray = (/pw&lat/)
res_pw@sfYArray = lev_pw
res_pw@cnFillOn = True ; turn on color fill
res_pw@gsnSpreadColors = True ; use full range of colo
res_pw@lbLabelAutoStride = True ; optimal labels
res_pw@trYReverse = True ; reverses y axis
;res_pw@trYLog = True

;second figure
;---------------------------
res_uw = True
res_uw@gsnDraw = False
res_uw@gsnFrame = False
res_uw@sfXArray = (/uwind&lat/)
res_uw@sfYArray = lev_uw
res_uw@trYReverse = True ; reverses y axis
res_uw@cnFillOn = True ; turn on color fill
res_uw@gsnSpreadColors = True ; use full range of colo
res_uw@lbLabelAutoStride = True ; optimal labels
;res_uw@trLog = True

;plot
;----
print("ploting...The CCA vectors must be reshaped")

plot_name=exp_name+"cca"
print("plot :"+plot_name)

wks=gsn_open_wks("eps",plot_name)
gsn_define_colormap(wks,"BlWhRe")
 
  ;
rcca_pw=onedtond(cca_pw(1,:),(/nlev_pw,nlat_pw/))
rcca_uw=onedtond(cca_uw(1,:),(/nlev_uw,nlat_uw/))
 
  
plot(0)=gsn_contour(wks,rcca_pw(:,:),res_pw)
plot(1)=gsn_contour(wks,rcca_uw(:,:),res_uw)

;create panel
;------------
gsn_panel(wks,plot,(/1,2/),False)

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Oct 22 12:16:02 2009

This archive was generated by hypermail 2.1.8 : Mon Oct 26 2009 - 15:01:33 MDT