;************************************************ ; bootstrap_correl_2.ncl ; ; Concepts illustrated: ; - Read SOI (Southern Oscilliation Index) Signal ; - Read gridded surface temperature ; - Calculate the monthly climatology and then anomalies ; - Calculate conventional cross correlation (SOI vs Anomalies) ; - Calculate 2.5% and 97.5% bounds via conventional statistical methods ; - Calculate bootstrapped cross correlations ; - Calculate 2.5% and 97.5% bounds via bootstrap methods ; - Create panel plot ;************************************************* 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 "./bootstrap.ncl_beta_640" ;************************************************* ;--- Read SOI data file an extract LSAT ;************************************************* dSOI = "/bcstor04/amita/" fSOI = "neind_sv_jjas_1979-2007.bin" ;timesereis pSOI = dSOI+fSOI UNDEF = -999.0 ntim = 29 ; UNDEF soi = new ( (/ntim/), float, "No_FillValue") soi = fbindirread(pSOI, 0, (/ntim/), "float") nsoi = dimsizes(soi) printVarSummary(soi) ymStrt = 1979 ymLast = 2007 print("ymStrt="+ymStrt+" ymLast="+ymLast) ;*********************************************************** ;--- Calculate some basic statistics for the SOI ; These statistics are not used .... just background statistics ;*********************************************************** soiStat4 = dim_stat4_n(soi, 0) ; 1st 4 moments of original sample ; for clarity; explicitly extract soiAvg = soiStat4(0) ; original sample mean soiStd = sqrt(soiStat4(1)) ; " sample std dev soiSkew = soiStat4(2) ; skewness; departure from symmetry soiKurt = soiStat4(3) ; kurtosis; relative to a normal distribution soiStdErr = soiStd/nsoi soiMed = dim_median_n(soi,0) ; median of original sample df = nsoi-1 p = 0.975 ; match default bootstrap 'p' (0.025, 0.975) tsoi = cdft_t(p,df) ; 2-sided soiLow = soiAvg-tsoi*soiStd ; normal low: 2.5% soiHi = soiAvg+tsoi*soiStd ; normal hi ; 97.5% ;*********************************************************** ;--- Read variable ; Calculate monthly climatology; Calculate monthly anomalies ;*********************************************************** diri = "/bcstor04/amita/" fili = "swe_mon_jan1979tomay2007smmrssmi.nc" pthi = diri+fili f = addfile(pthi, "r") y = f->swe printVarSummary(y) yAnom = new ( (/29, 60, 360/), float, UNDEF) do i=0,28 yAnom(i,:,:) = y(i*12,:,:) end do copy_VarMeta(y(0:28,:,:),yAnom) printVarSummary(yAnom) dimy = dimsizes(yAnom) ntim = dimy(0) if (ntim.ne.nsoi) then print("time mismatch: ntim="+ntim+" nsoi="+nsoi) exit end if ;*********************************************************** ;--- Compute conventional sample correlation('r') at each grid point ; Assumes bivariate normal distribution ; 'r' distribution is skewed, hence, use Fischer z-transform ;*********************************************************** r = escorc_n(soi,yAnom,0,0) r@long_name = "correlation: escorc" df = nsoi-2 rse = sqrt(1-r^2)/sqrt(df) rse@long_name = "correlation: escorc: std. error" z = 0.5*log((1+r)/(1-r)) z@long_name = "Fischer z-transform" zse = 1.0/sqrt(nsoi-3) zse@long_name = "standard error of z statistic" ;*********************************************************** ;--- Add coordinate meta data ;*********************************************************** copy_VarCoords(y(0,:,:), r) ; create (copy) spatial coordinates printVarSummary(r) copy_VarCoords(y(0,:,:), rse) ; create (copy) spatial coordinates printVarSummary(rse) copy_VarCoords(y(0,:,:), z) ; create (copy) spatial coordinates printVarSummary(z) ;*********************************************************** ;--- Conventional (normal) estimates; use z-transformed values ;--- Get t-value for specified 'ptst' : Add meta data ;*********************************************************** ptst = 0.975 ; match default bootstrap 'p' (0.025, 0.975) zval = cdft_t(ptst,9999) ; 2-sided; use 'big' number for normal dist zLow = z-zval*zse ; z-normal low: 2.5% zHi = z+zval*zse ; z-normal hi ; 97.5% ; transform z back to 'r' space rLow = tanh(zLow) ; =(exp(2*zlow)-1)/(exp(2*zlow)+1) rHi = tanh(zHi) ; =(exp(2*zhi )-1)/(exp(2*zhi )+1) rLow@long_name = "Correlation: 2.5%" rHi@long_name = "Correlation: 97.5%" tval = r*sqrt(df/(1-r^2)) ; conventional t-value psoi = student_t(tval,df) ; probability psoi@long_name = "Conventional Probability" copy_VarCoords(y(0,:,:), rLow); create (copy) spatial coordinates copy_VarCoords(y(0,:,:), rHi ) copy_VarCoords(y(0,:,:), psoi) printVarSummary(rLow) printVarSummary(rHi) printVarSummary(psoi) ;*********************************************************** ;--- BOOTSTRAP ; Extract from returned 'list' vatiable; Add meta data ;*********************************************************** N = ntim ; convenience n = N ; no subsampling nBoot = 500 opt = False if (n.ne.N) then opt = True opt@sample_size = n end if Boot_r = bootstrap_correl(soi, yAnom, nBoot, (/0,0/), opt) rBoot = Boot_r[0] ; bootstrapped correlations rBootAvg = Boot_r[1] ; average of bootstrapped correlations rBootStd = Boot_r[2] ; std dev of bootstrapped correlations delete(Boot_r) ; no longer needed rBootLow = bootstrap_estimate(rBoot, 0.025, False) ; 2.5% lower confidence bound rBootMed = bootstrap_estimate(rBoot, 0.500, False) ; 50.0% median of bootstrapped estimates rBootHi = bootstrap_estimate(rBoot, 0.975, False) ; 97.5% upper confidence bound rBootLow@long_name = "bootstrap r: "+rBootLow@estimate rBootMed@long_name = "bootstrap r: Median" rBootHi@long_name = "bootstrap r: "+rBootHi@estimate printVarSummary(rBoot) printMinMax(rBoot, 0) printVarSummary(rBootAvg) printMinMax(rBootAvg, 0) ;+++++++++++++++++++++++++++++++++++ ; PLOTS ;+++++++++++++++++++++++++++++++++++ pltDir = "./" ;pltName = "BootCor_SOI_"+N+"_"+n+"_"+nBoot pltName = "bootstrap_correl" pltPath = pltDir+pltName pltType = "x11" wks = gsn_open_wks (pltType,pltPath) ;*************************************************************** ;--- create histogram for the SOI_SIGNAL ;*************************************************************** resh = True resh@gsnDraw = False resh@gsnFrame = False resh@gsnHistogramNumberOfBins = 21 resh@gsFillColor = "green" ; resh@tiMainString = "Original: N="+N hstSoi = gsn_histogram(wks, soi ,resh) ;*************************************************************** ;--- text object SOI_SIGNAL ;*************************************************************** txres = True txres@txFont = "helvetica-bold" txres@txFontHeightF = 0.0150 textSoi = (/" Mean="+sprintf("%5.1f", soiAvg) +"~C~"+ \ " Std="+sprintf("%5.1f", soiStd) +"~C~"+ \ "StdErr="+sprintf("%5.1f", soiStdErr) +"~C~"+ \ " Skew="+sprintf("%5.1f", soiSkew) +"~C~"+ \ " Kurt="+sprintf("%5.1f", soiKurt) +"~C~"+ \ " Low="+sprintf("%5.1f", soiLow) +"~C~"+ \ "Median="+sprintf("%5.1f", soiMed) +"~C~"+ \ " Hi="+sprintf("%5.1f", soiHi ) /) txBoxSoi = gsn_create_text(wks,textSoi,txres) amres = True amres@amParallelPosF = 0.35 ; move legend to the right amres@amOrthogonalPosF = -0.275 ; move the legend up annoSoi = gsn_add_annotation(hstSoi, txBoxSoi, amres) ; Attach string to plot draw(hstSoi) frame(wks) ;*************************************************************** ;--- Global correlations of SOI_SIGNAL and monthly anomalies ;*************************************************************** rescn = True rescn@gsnDraw = False rescn@gsnFrame = False rescn@mpFillOn = False ; no need rescn@cnLevelSelectionMode= "ManualLevels" ; manual set levels rescn@cnMinLevelValF = -1.0 rescn@cnMaxLevelValF = 1.0 rescn@cnLevelSpacingF = 0.2 rescn@cnFillOn = True ; color fill plot rescn@cnLinesOn = False rescn@cnLineLabelsOn = False rescn@cnInfoLabelOn = False rescn@cnFillMode = "RasterFill" rescn@gsnSpreadColors = True rescn@lbLabelBarOn = False ; turn off individual label bars rescn@gsnStringFontHeightF= 0.02 rescn@mpMinLatF = 30 ; specify min lat rescn@gsnPolar = "NH" ; specify the hemisphere ;---Create arrays to hold series of plots plots = new(6,graphic) gsn_define_colormap(wks,"MPL_BrBG") plots(0) = gsn_csm_contour_map_polar(wks,rLow,rescn) plots(1) = gsn_csm_contour_map_polar(wks,rBootLow,rescn) plots(2) = gsn_csm_contour_map_polar(wks,r ,rescn) plots(3) = gsn_csm_contour_map_polar(wks,rBootMed,rescn) plots(4) = gsn_csm_contour_map_polar(wks,rHi ,rescn) plots(5) = gsn_csm_contour_map_polar(wks,rBootHi ,rescn) ;---Resources for paneling pres = True ; modify the panel plot pres@txString = "SOI-TempAnom Cor: "+ymStrt+"-"+ymLast+ \ ": Conventional/BootStrap: nBoot="+nBoot+": N="+N+", n="+n pres@txFontHeightF = 0.0175 ; default 0.05 pres@gsnPanelLabelBar = True ; add common colorbar gsn_panel(wks,plots,(/3,2/),pres)