[ncl-talk] ttest: probabilities with very small values
Sri.durgesh Nandini-Weiss
sri.durgesh.nandini-weiss at uni-hamburg.de
Thu Nov 14 07:45:57 MST 2019
Hello
Performing ttest on the 2 original samples yields the same result, with
no warning nor error.
The original samples mean:
(0) min=2.11077 max=2.11077
(0) min=0.959119 max=0.959119
and now the p value is (0) min=2.29254e-10 max=2.29254e-10
BUT: the next step is to generate PDFs of each sample distribution
(normally this is always done on anomalies), and the PDFs are simple
strange.
I have attached my script and the pdf plot here when using the (1)
original sample+ pdf and when using the (2) anomalies+pdf.
As you see, the one using anomaly data is the correct way for
calculating and plotting the pdf, but i cannot understand how this will
be archived if the means are so small.
Would be grateful for any help.
Sri
On 11/14/19 2:43 PM, Dennis Shea wrote:
> I have not looked at the script carefully. Some comments
>
> [1] NCL's *ttest*
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml> is
> essentially the same as that in *Numerical Recipes*
> <https://websites.pmc.ucsc.edu/~fnimmo/eart290c_17/NumericalRecipesinF77.pdf>
> [*pg 610; Chapter 14*].
> [2] re: warning:*ttest:* encountered 6 cases where s1 and/or s2 were
> less than 2.
> s1 and s2 refer to sample sizes*
> *
> [3] *equiv_sample_size*
> <https://www.ncl.ucar.edu/Document/Functions/Built-in/equiv_sample_size.shtml>
> **Apparently, this is returning very small sample sizes*?*due to very
> high [near 1.0] autocorrelation*?
> *
> **For 'fun', try running the script**using the full sample size*
> *
>
> On Wed, Nov 13, 2019 at 3:51 AM Sri.durgesh Nandini-Weiss via ncl-talk
> <ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>> wrote:
>
> Good morning everyone,
>
> I was wondering to get some advice on the below.
>
> I have a working script for computing ttest, which is computed on
> anomalies, however, the values are incredibly small-
>
> at the significance of 0.05 returned by ttest would yield 95% for
> alpha.
>
> Probablities : min=0 max=0.000435114
> SO im wondering if it actually makes sense to try to plot it. I have
> only one warning:
>
> warning:ttest: encountered 6 cases where s1 and/or s2 were less
> than 2.
> Output set to missing values in these cases.
>
> Heres my script:
>
> ;========================================================================================
> ; Concepts illustrated in my script:
> sri.durgesh.nandini-weiss at uni-hamburg.de
> <mailto:sri.durgesh.nandini-weiss at uni-hamburg.de>
> ; The attached calculates the ensemble (N=100) mean wave height,
> standard deviation, skewness and kurtosis between 1981-2006 and
> 2080-2100
> ; Then it plots ensembl mean higher moments in statistics by
> computing
> 20yrs anomalies 4 (3) moments in panel:
> ; (a) Monthly climatology and monthly anomalies
> ; (b) Ensemble spread (std. deviation)
> ; (c) skewness and (d) kurtosis
>
> ;==============================================================
> ; Open the file: Read only the user specified period
> ; =============================================================
>
> f = addfile ("anom_hist_zo.nc <http://anom_hist_zo.nc>", "r")
> hist_anom = f->hist_anom
> printVarSummary(hist_anom) ;[time | 240] x [ens | 100] x [lat | 45] x
> [lon | 90]
>
> f1 = addfile ("anom_rcp45_zo.nc
> <http://anom_rcp45_zo.nc>", "r")
> rcp45_anom = f1->rcp45_anom
> printVarSummary(rcp45_anom) ;[time | 240] x [ens | 100] x [lat
> | 45]
> x [lon | 90]
>
>
> dimx = dimsizes(hist_anom)
> ntim = dimx(0) ; 240
> nens = dimx(1) ; 100
> nlat = dimx(2) ; 45
> mlon = dimx(3) ; 90
>
> nmos = 12
> nyrs = ntim/nmos ; 20
> printVarSummary(nyrs)
>
> ;==================================================================
> ; Compute first four moments (average, variance, skewness, and
> kurtosis)
> ; over the time and ensemble dimension together: for Hist, and
> then RCP45
> ;==================================================================
>
> aveX = new((/nlat,mlon/), typeof(hist_anom),
> hist_anom at _FillValue)
> varX = new((/nlat,mlon/), typeof(hist_anom),
> hist_anom at _FillValue)
> xSkwMonEns = new((/nlat,mlon/), typeof(hist_anom),
> hist_anom at _FillValue)
> xKurMonEns = new((/nlat,mlon/), typeof(hist_anom),
> hist_anom at _FillValue)
>
> do nmo=0,nmos-1
> work := reshape(hist_anom(nmo::nmos,:,:,:)
> ,(/nyrs*nens,nlat,mlon/))
> if (nmo.eq.0) then
> printVarSummary(work) ;
> (2000,45,90)
> end if
> xStat = dim_stat4_n(work,0 )
> end do
> printVarSummary(xStat) ;
> (4,nlat,mlon)
> aveX = xStat(0,:,:) ; (nlat,
> mlon)
> varX = xStat(1,:,:)
> xSkwMonEns = xStat(2,:,:)
> xKurMonEns = xStat(3,:,:)
>
> copy_VarCoords(hist_anom(0,0,:,:), aveX)
> aveX at long_name = "Monthly Climatology (Average)"
> aveX at LONG_NAME = "Monthly Climatology over all Ensemble Members"
>
> copy_VarCoords(hist_anom(0,0,:,:), varX)
> varX at long_name = "Monthly Interannual Variability (sample
> variance)"
> varX at LONG_NAME = "Monthly Interannual Variability over all
> Ensemble
> Members"
>
> copy_VarCoords(hist_anom(0,0,:,:), xSkwMonEns)
> xSkwMonEns at long_name = "Monthly Skewness"
> copy_VarCoords(hist_anom(0,0,:,:), xKurMonEns)
> xKurMonEns at long_name = "Monthly Kurtosis"
>
> print("============")
>
> ;================================================================================
> ; dim average on month dimension; does this change the plotted
> results
> and prob?
> ;================================================================================
>
> printVarSummary(aveX) ; (lat,lon)
> printMinMax(aveX,0)
> printVarSummary(varX) ; (lat,lon)
> printMinMax(varX,0)
> printVarSummary(xSkwMonEns) ; (lat,lon)
> printMinMax(xSkwMonEns,0)
> printVarSummary(xKurMonEns) ; (lat,lon)
> printMinMax(xKurMonEns,0)
> print("============")
>
> ;==================================================================
> ; Now repeat for RCP45
> ;==================================================================
>
>
> aveY = new((/nlat,mlon/), typeof(rcp45_anom),
> rcp45_anom at _FillValue)
> varY = new((/nlat,mlon/), typeof(rcp45_anom),
> rcp45_anom at _FillValue)
> xSkwMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom),
> rcp45_anom at _FillValue)
> xKurMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom),
> rcp45_anom at _FillValue)
>
> do nmo=0,nmos-1
> work2 := reshape(rcp45_anom(nmo::nmos,:,:,:)
> ,(/nyrs*nens,nlat,mlon/))
> if (nmo.eq.0) then
> printVarSummary(work2) ;
> (2000,45,90)
> end if
> xStat2 = dim_stat4_n(work2,0 )
> end do
> printVarSummary(xStat2) ;
> (4,nlat,mlon)
> aveY = xStat2(0,:,:) ;
> (nlat, mlon)
> varY = xStat2(1,:,:)
> xSkwMonEns2 = xStat2(2,:,:)
> xKurMonEns2 = xStat2(3,:,:)
>
> copy_VarCoords(rcp45_anom(0,0,:,:), aveY)
> aveY at long_name = "Monthly Climatology (Average)"
> aveY at LONG_NAME = "Monthly Climatology over all Ensemble Members"
>
> copy_VarCoords(rcp45_anom(0,0,:,:), varY)
> varY at long_name = "Monthly Interannual Variability (sample
> variance)"
> varY at LONG_NAME = "Monthly Interannual Variability over all
> Ensemble
> Members"
>
> copy_VarCoords(rcp45_anom(0,0,:,:), xSkwMonEns2)
> xSkwMonEns2 at long_name = "Monthly Skewness"
> copy_VarCoords(rcp45_anom(0,0,:,:), xKurMonEns2)
> xKurMonEns2 at long_name = "Monthly Kurtosis"
>
> print("============")
>
> printVarSummary(aveY) ; (lat,lon)
> printMinMax(aveY,0)
> printVarSummary(varY) ; (lat,lon)
> printMinMax(varY,0)
> printVarSummary(xSkwMonEns2) ; (lat,lon)
> printMinMax(xSkwMonEns2,0)
> printVarSummary(xKurMonEns2) ; (lat,lon)
> printMinMax(xKurMonEns2,0)
> print("============")
>
> ;================================================================================
> ; equiv_sample_size: Estimates the number of independent values of a
> series of correlated observations.
> ; Specify a critical significance level to test the lag-one
> auto-correlation coefficient.
> ; monthly climatolgy perharps not enough data points to calculate
> ;=================================================================================
> ; sigr = 0.01 or 0.05 If prob < sigr, then the null hypothesis
> (means
> are from the same population)
> ; is rejected and the alternative hypothesis is accepted.
> ;================================================================================
>
> sigr = 0.05 ;all areas
> whose
> value > 95
>
> sX = new(dimsizes(aveX),typeof(aveX),aveX at _FillValue)
> printVarSummary(sX)
>
> sX = equiv_sample_size (hist_anom(lat|:,lon|:,ens|:0,time|:),
> sigr,0) ; (time,lat,lon)
> copy_VarMeta(aveX,sX)
> printMinMax(sX,0)
> printVarSummary(sX) ; (lat,lon)
> ;print(sX)
>
> sY = new(dimsizes(aveX),typeof(aveX),aveX at _FillValue)
> printVarSummary(sY)
>
> sY = equiv_sample_size (rcp45_anom(lat|:,lon|:,ens|:0,time|:),
> sigr,0) ; (time,lat,lon)
> copy_VarMeta(aveY,sY)
> printMinMax(sY,0)
> printVarSummary(sY) ; (lat,lon)
> ;print(sY)
>
>
> ;==================================================================================
> ; Use "ttest" to compute the probabilities.The returned
> probability is
> two-tailed.
> ;==================================================================================
>
> iflag= True ; population
> variance similar?
> prob = new(dimsizes(varY),typeof(varY),varY at _FillValue)
> printVarSummary(prob)
>
> ;==================================================================================
> ; A significance of 0.05 returned by ttest would yield 95% for alpha.
> ;==================================================================================
>
> prob = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, iflag, False))
> copy_VarCoords(aveY, prob)
> prob at long_name = " Probablities"
> printVarSummary(prob) ; (lat,lon)
> printMinMax(prob,0) ; min=0 max=100
> ;print(prob)
>
> ;==================================================================================
>
> Best
>
> Sri
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
--
Dr. Sri Nandini-Weiß
Research Scientist
Universität Hamburg
Center for Earth System Research and Sustainability (CEN)
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
Bundesstrasse 53, 20146 Hamburg
Tel: +49 (0) 40 42838 7472
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0001.html>
-------------- next part --------------
;========================================================================================
; Concepts illustrated in my script: sri.durgesh.nandini-weiss at uni-hamburg.de
; The attached calculates the Basic Extreme Value Statistics to analyse extreme deviations from the median of probability distributions for 1981-2006 and 2080-2100
; Then it plots ensembl mean higher moments in statistics by remove a climatological seasonal cycle: as t-test (mean) and f-test (variance):
; (i)Calculates the probability (PDF) and cumulative (CDF) distribution functions of the Generalized Extreme Value (GEV) distribution given the shape, scale and location parameters.
; (a) Ensemble mean and (b) spread (variance)
; (c) skewness and (d) kurtosis (e) 5th and 95th percentiles of the respective dsitributions.
; calculate and plot PDFs of simulated sea surface height anomalies (cm) for island regions in the tropical Pacific.
; but first the distribution functions are presented for
; (a) the globe
; (b) northern hermisphere extra tropics (30N-60N)
; (c) tropics (30S-30N)
; (d) souther hermipere extra tropics (30-90S)
; Aim: To understand the shifting probability distribution of SSH under rcp45 and rcp85.
; For this the seasonal cycles of the respective periods should be removed otherwise probably not Gaussian.
; One cycle per period not per member.
; extreme value analyses focus on high extremes: maxima (SSH)
;========================================================================================
;==============================================================
; Open the file: Read only the user specified period
; =============================================================
begin
yrStrt = 1986
yrLast = 2006
f = addfile("zo_1850-2005_ens_1-100.nc","r") ; open netcdf file
time = f->time
printVarSummary(time)
time = time-31
YYYY = cd_calendar(time,-1)/100 ; zo:[time|1872] x [ens|100] x [depth|1] x [lat|45] x [lon|90]
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
x = f->zo(iYYYY,:,0,:,:)
x at _FillValue = -9.96921e+36
hist_anom = x
printVarSummary(hist_anom)
yrStrt4 = 2080
yrLast4 = 2100
f4 = addfile("zo_RCP45_2006-2100_ens_1-100.nc","r") ; open netcdf file
time4 = f4->time
time4 = time4-31
printVarSummary(time4)
YYYY4 = cd_calendar(time4,-1)/100 ; zo:[time|1872] x [ens|100] x [depth|1] x [lat|45] x [lon|90]
iYYYY4 = ind(YYYY4.ge.yrStrt4 .and. YYYY4.le.yrLast4)
x2 = f4->zo(iYYYY4,:,0,:,:)
printVarSummary(x2)
x at _FillValue = -9.96921e+36
rcp45_anom = x2
printVarSummary(rcp45_anom)
dimx = dimsizes(hist_anom)
ntim = dimx(0) ; 240
nens = dimx(1) ; 100
nlat = dimx(2) ; 45
nlon = dimx(3) ; 90
lat1=nlat
lat2=nlat
nmos = 12
nyrs = ntim/nmos ; 20
printVarSummary(nyrs)
; =============================================================
regions=(/"Global","30N-90N","30S-30N","90S-30S","10-20"/)
do ireg=0,dimsizes(regions)-1
region=regions(ireg)
print(ireg+": "+region)
if(region.eq."global") then
lat1=-90
lat2=90
maxY=0.08
precisionY=1
end if
if(region.eq."30N-90N") then
lat1=30
lat2=90
maxY=0.07
precisionY=1
end if
if(region.eq."30S-30N") then
lat1=-30
lat2=30
maxY=0.13
precisionY=2
end if
if(region.eq."90S-30S") then
lat1=-90
lat2=-30
maxY=0.07
precisionY=1
end if
if(region.eq."20S-20N") then
lat1=-20
lat2=20
maxY=0.12
precisionY=2
end if
if (ireg.gt.0) then
delete(anom1_nomiss)
delete(anom1_permsort)
delete(anom2_nomiss)
delete(anom2_permsort)
delete( hist_anom1)
delete( rcp45_anom1)
delete(anom1_1d)
delete(anom2_1d)
delete(res)
end if
;==================================================================
print ("Data read in - start calculations")
;==================================================================
; Region domain lat and lon
;==================================================================
hist_anom1=hist_anom(:,:,{lat1:lat2},:) ;(ntim, mlon)
printVarSummary(hist_anom1)
rcp45_anom1=rcp45_anom(:,:,{lat1:lat2},:) ;(ntim, mlon)
printVarSummary(rcp45_anom1)
;=================================================================
anom1_1d = ndtooned(hist_anom1) ;to carry out the dim_stat4_n on one dimension.
printVarSummary(anom1_1d)
anom2_1d = ndtooned(rcp45_anom1)
printVarSummary(anom2_1d)
;==================================================================
TTstat1 = dim_stat4(anom1_1d)
printVarSummary(TTstat1) ; (4,?)
TTstat2 = dim_stat4(anom2_1d)
printVarSummary(TTstat2) ; (4,?)
print("============")
;p=addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
;lsmask=landsea_mask(p->LSMASK,ts3_diff&lat,ts3_diff&lon)
;diff1_lnd=ts3_diff
;diff1_lnd=mask(ts3_diff,lsmask.le.0.0,True)
;delete(lsmask)
;ttest
sigr = 0.05 ;all areas whose value > 95
aveX = avg (anom1_1d)
printMinMax(aveX,0)
aveY = avg (anom2_1d)
printMinMax(aveY,0)
varX = variance (anom1_1d)
varY = variance (anom2_1d)
;==================================================================
; if the sample means in two groups are identical, the p-values of a t-test is 1.
; degrees of freedom as number of independent values in series
;==================================================================
anom1_1d(ind(ismissing(anom1_1d))) = -9999 ;(too many missing values - remove)
anom2_1d(ind(ismissing(anom2_1d))) = -9999
sX = equiv_sample_size (anom1_1d(ind(anom1_1d.ne.-9999)), 0.10,0)
sY = equiv_sample_size (anom2_1d(ind(anom2_1d.ne.-9999)), 0.10,0)
printVarSummary(sY)
iflag= True
;t_prob2 = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, iflag, False) )
t_prob2 = ttest(aveX,varX,sX, aveY,varY,sY, iflag, False)
printVarSummary(t_prob2)
printMinMax(t_prob2,0)
;f_prob = 100.*(1. -ftest(varX,sX,varY,sY,0)) ; test significance of variances!!
f_prob = ftest(varX,sX,varY,sY,0)
printVarSummary(f_prob)
printMinMax(f_prob,0)
;==================================================================
print(sX+" "+sY+" "+t_prob2+" "+f_prob)
; percentiles (requires sorting the array)
anom1_nomiss=anom1_1d(ind(anom1_1d.ne.-9999))
anom1_permsort = dim_pqsort(anom1_nomiss,2) ;; sorted array is anom1_nomiss!!
anom1_p05=anom1_nomiss(floattoint(0.05*dimsizes(anom1_nomiss))) ;(rank_p5)
anom1_p10=anom1_nomiss(floattoint(0.1*dimsizes(anom1_nomiss)))
anom1_p90=anom1_nomiss(floattoint(0.9*dimsizes(anom1_nomiss)))
anom1_p95=anom1_nomiss(floattoint(0.95*dimsizes(anom1_nomiss)))
printVarSummary(anom1_p05)
anom2_nomiss=anom2_1d(ind(anom2_1d.ne.-9999))
anom2_permsort = dim_pqsort(anom2_nomiss,2) ; sorted array is anom1_nomiss!!
anom2_p05=anom2_nomiss(floattoint(0.05*dimsizes(anom2_nomiss))) ;(rank_p5)
anom2_p10=anom2_nomiss(floattoint(0.1*dimsizes(anom2_nomiss)))
anom2_p90=anom2_nomiss(floattoint(0.9*dimsizes(anom2_nomiss)))
anom2_p95=anom2_nomiss(floattoint(0.95*dimsizes(anom2_nomiss)))
print(anom2_p05+" "+anom2_p10+" "+anom2_p90+" "+anom2_p95)
print("calculate PDFs=========================")
;================================================================================
; calculate empirical PDF
;================================================================================
opt = True
opt at bin_min = -25.
opt at bin_max = 25.
;opt at bin_nice = True
pdf1 = pdfx(hist_anom1, 100, opt)
printVarSummary(pdf1)
pdf2 = pdfx(rcp45_anom1, 100, opt)
printVarSummary(pdf2)
;================================================================================
; calculate y-value for percentiles
;================================================================================
p05pos=ind(pdf1 at bin_center.ge.anom1_p05)
y05_1 = (anom1_p05-pdf1 at bin_center(p05pos(0)-1)) / (pdf1 at bin_center(p05pos(1)-1)-pdf1 at bin_center(p05pos(0)-1)) * (pdf1(p05pos(1)-1)-pdf1(p05pos(0)-1)) + pdf1(p05pos(0)-1)
p95pos=ind(pdf1 at bin_center.ge.anom1_p95)
y95_1 = (anom1_p95-pdf1 at bin_center(p95pos(0)-1)) / (pdf1 at bin_center(p95pos(1)-1)-pdf1 at bin_center(p95pos(0)-1)) * (pdf1(p95pos(1)-1)-pdf1(p95pos(0)-1)) + pdf1(p95pos(0)-1)
delete(p05pos)
delete(p95pos)
p05pos=ind(pdf2 at bin_center.ge.anom2_p05)
y05_2 = (anom2_p05-pdf2 at bin_center(p05pos(0)-1)) / (pdf2 at bin_center(p05pos(1)-1)-pdf2 at bin_center(p05pos(0)-1)) * (pdf2(p05pos(1)-1)-pdf2(p05pos(0)-1)) + pdf2(p05pos(0)-1)
p95pos=ind(pdf2 at bin_center.ge.anom2_p95)
y95_2 = (anom2_p95-pdf2 at bin_center(p95pos(0)-1)) / (pdf2 at bin_center(p95pos(1)-1)-pdf2 at bin_center(p95pos(0)-1)) * (pdf2(p95pos(1)-1)-pdf2(p95pos(0)-1)) + pdf2(p95pos(0)-1)
delete(p05pos)
delete(p95pos)
;================================================================================
; prepare data for plotting
;================================================================================
nVar = 2
nBin = pdf1 at nbins ; retrieve the number of bins
xx = new ( (/nVar, nBin/), typeof(pdf1))
xx(0,:) = pdf1 at bin_center ; assign appropriate "x" axis values
xx(1,:) = pdf2 at bin_center
yy = new ( (/nVar, nBin/), typeof(pdf1))
yy(0,:) = (/ pdf1 /) * 0.01 ;; convert % to absolut
yy(1,:) = (/ pdf2 /) * 0.01
;================================================================================
; Plot PDFs...seperate plots
;================================================================================
plotfile="Hist_RCP45_"+region+"_PDF"
wks_type = "pdf"
wks = gsn_open_wks(wks_type,plotfile)
res = True
res at gsnMaximize =True
res at vpHeightF = 0.7 ; change aspect ratio of plot
res at vpWidthF = 1.
res at tmYLAutoPrecision = False ; no auto precision
res at tmYLPrecision = 2 ; set the precision
res at xyLineThicknesses = (/5.,5./)
res at xyDashPattern = 0 ; dashed/solid(=0)
;res at trXMinF = -15 ; set minimum X-axis value
;res at trXMaxF = 15 ; set maximum X-axis value
;res at trYMaxF = 0.13 ; set maximum Y-axis value
res at tiXAxisFontHeightF = 0.02
res at tiYAxisFontHeightF = 0.02
res at tiYAxisString = "probability %" ; axis string
res at tiXAxisString = "SSH anomaly [cm]" ; axis string
res at xyLineColors = (/"blue","red"/)
res at gsnXRefLineDashPattern =1
res at gsnXRefLine = 0.0
res at gsnDraw = False ; don't draw so text can be added
res at gsnFrame = False
res at gsnStringFontHeightF = 0.02
res at gsnCenterString = "Region:" +region
plot = gsn_csm_xy (wks, xx, yy, res)
;================================================================================
; vertical lines for percentiles
;================================================================================
plines = new(4,graphic)
plres=True
plres at gsLineColor = "blue"
plres at gsLineThicknessF = 5.
plines(0) = gsn_add_polyline(wks, plot, (/anom1_p05,anom1_p05/), (/0,0.01*y05_1/), plres)
plines(1) = gsn_add_polyline(wks, plot, (/anom1_p95,anom1_p95/), (/0,0.01*y95_1/), plres)
plres at gsLineColor = "red"
plines(2) = gsn_add_polyline(wks, plot, (/anom2_p05,anom2_p05/), (/0,0.01*y05_2/), plres)
plines(3) = gsn_add_polyline(wks, plot, (/anom2_p95,anom2_p95/), (/0,0.01*y95_2/), plres)
;================================================================================
; add text
;================================================================================
txres = True
;txres at txFont = "duplex_roman"
txres at txJust = "CenterLeft"
txres at txFontHeightF = 0.022
;txres at txFuncCode = "~"
txres at txFontColor = "black"
;; with axis lables
gsn_text_ndc (wks,"Significance", 0.19, 0.78,txres)
gsn_text_ndc (wks,"T-Test: "+decimalPlaces(t_prob2,3,True), 0.19, 0.745,txres)
gsn_text_ndc (wks,"F-Test: "+decimalPlaces(f_prob,3,True), 0.191, 0.71,txres)
txres at txFontColor = "blue"
gsn_text_ndc(wks,"1986-2006",0.19,0.645,txres)
gsn_text_ndc (wks,"Mean = "+decimalPlaces(TTstat1(0),3,True), 0.192, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat1(1),3,True), 0.192, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat1(2),3,True), 0.192, 0.54 ,txres)
;gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat1(3),3,True), 0.192, 0.505 ,txres)
txres at txFontColor = "red"
gsn_text_ndc(wks,"RCP4.5:1981-2100",0.698,0.645,txres)
gsn_text_ndc (wks,"Mean = "+decimalPlaces(TTstat2(0),3,True), 0.70, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat2(1),3,True), 0.70, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat2(2),3,True), 0.70, 0.54 ,txres)
;gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat2(3),3,True), 0.70, 0.505 ,txres)
draw (plot)
frame (wks)
end do ; loop regions
end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_10-20_PDF.pdf
Type: application/pdf
Size: 95850 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0007.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_30N-90N_PDF.pdf
Type: application/pdf
Size: 93475 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0008.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_90S-30S_PDF.pdf
Type: application/pdf
Size: 98491 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0009.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_Global_PDF.pdf
Type: application/pdf
Size: 90506 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0010.pdf>
-------------- next part --------------
;========================================================================================
; Concepts illustrated in my script: sri.durgesh.nandini-weiss at uni-hamburg.de
; The attached calculates the ensemble (N=100) mean wave height, standard deviation, skewness and kurtosis between 1981-2006 and 2080-2100
; Then it plots ensembl mean higher moments in statistics by computing 20yrs anomalies 4 (3) moments in panel:
; (a) Monthly climatology and monthly anomalies
; (b) Ensemble spread (std. deviation)
; (c) skewness and (d) kurtosis
;==============================================================
; Open the file: Read only the user specified period
; =============================================================
f = addfile ("anom_hist_zo.nc", "r")
hist_anom = f->hist_anom
printVarSummary(hist_anom)
f1 = addfile ("anom_rcp45_zo.nc", "r")
rcp45_anom = f1->rcp45_anom
printVarSummary(rcp45_anom)
dimx = dimsizes(hist_anom)
ntim = dimx(0) ; 240
nens = dimx(1) ; 100
nlat = dimx(2) ; 45
mlon = dimx(3) ; 90
nmos = 12
nyrs = ntim/nmos ; 20
printVarSummary(nyrs)
;==================================================================
; Compute first four moments (average, variance, skewness, and kurtosis)
; over the time and ensemble dimension together: for Hist, and then RCP45
;==================================================================
aveX = new((/nlat,mlon/), typeof(hist_anom), hist_anom at _FillValue)
varX = new((/nlat,mlon/), typeof(hist_anom), hist_anom at _FillValue)
xSkwMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom at _FillValue)
xKurMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom at _FillValue)
do nmo=0,nmos-1
work := reshape(hist_anom(nmo::nmos,:,:,:) ,(/nyrs*nens,nlat,mlon/))
if (nmo.eq.0) then
printVarSummary(work) ; (2000,45,90)
end if
xStat = dim_stat4_n(work,0 )
end do
printVarSummary(xStat) ; (4,nlat,mlon)
aveX = xStat(0,:,:) ; (nlat, mlon)
varX = xStat(1,:,:)
xSkwMonEns = xStat(2,:,:)
xKurMonEns = xStat(3,:,:)
copy_VarCoords(hist_anom(0,0,:,:), aveX)
aveX at long_name = "Monthly Climatology (Average)"
aveX at LONG_NAME = "Monthly Climatology over all Ensemble Members"
copy_VarCoords(hist_anom(0,0,:,:), varX)
varX at long_name = "Monthly Interannual Variability (sample variance)"
varX at LONG_NAME = "Monthly Interannual Variability over all Ensemble Members"
copy_VarCoords(hist_anom(0,0,:,:), xSkwMonEns)
xSkwMonEns at long_name = "Monthly Skewness"
copy_VarCoords(hist_anom(0,0,:,:), xKurMonEns)
xKurMonEns at long_name = "Monthly Kurtosis"
print("============")
;================================================================================
; dim average on month dimension; does this change the plotted results and prob?
;================================================================================
;aveX = dim_avg_n_Wrap(aveX,0)
;varX = dim_avg_n_Wrap(varX,0)
;xSkwMonEns = dim_avg_n_Wrap(xSkwMonEns,0)
;xKurMonEns = dim_avg_n_Wrap(xKurMonEns,0)
printVarSummary(aveX) ; (lat,lon)
printMinMax(aveX,0)
printVarSummary(varX) ; (lat,lon)
printMinMax(varX,0)
printVarSummary(xSkwMonEns) ; (lat,lon)
printMinMax(xSkwMonEns,0)
printVarSummary(xKurMonEns) ; (lat,lon)
printMinMax(xKurMonEns,0)
print("============")
;==================================================================
; Now repeat for RCP45
;==================================================================
aveY = new((/nlat,mlon/), typeof(rcp45_anom), rcp45_anom at _FillValue)
varY = new((/nlat,mlon/), typeof(rcp45_anom), rcp45_anom at _FillValue)
xSkwMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom), rcp45_anom at _FillValue)
xKurMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom), rcp45_anom at _FillValue)
do nmo=0,nmos-1
work2 := reshape(rcp45_anom(nmo::nmos,:,:,:) ,(/nyrs*nens,nlat,mlon/))
if (nmo.eq.0) then
printVarSummary(work2) ; (2000,45,90)
end if
xStat2 = dim_stat4_n(work2,0 )
end do
printVarSummary(xStat2) ; (4,nlat,mlon)
aveY = xStat2(0,:,:) ; (nlat, mlon)
varY = xStat2(1,:,:)
xSkwMonEns2 = xStat2(2,:,:)
xKurMonEns2 = xStat2(3,:,:)
copy_VarCoords(rcp45_anom(0,0,:,:), aveY)
aveY at long_name = "Monthly Climatology (Average)"
aveY at LONG_NAME = "Monthly Climatology over all Ensemble Members"
copy_VarCoords(rcp45_anom(0,0,:,:), varY)
varY at long_name = "Monthly Interannual Variability (sample variance)"
varY at LONG_NAME = "Monthly Interannual Variability over all Ensemble Members"
copy_VarCoords(rcp45_anom(0,0,:,:), xSkwMonEns2)
xSkwMonEns2 at long_name = "Monthly Skewness"
copy_VarCoords(rcp45_anom(0,0,:,:), xKurMonEns2)
xKurMonEns2 at long_name = "Monthly Kurtosis"
print("============")
printVarSummary(aveY) ; (lat,lon)
printMinMax(aveY,0)
printVarSummary(varY) ; (lat,lon)
printMinMax(varY,0)
printVarSummary(xSkwMonEns2) ; (lat,lon)
printMinMax(xSkwMonEns2,0)
printVarSummary(xKurMonEns2) ; (lat,lon)
printMinMax(xKurMonEns2,0)
print("============")
;================================================================================
; equiv_sample_size: Estimates the number of independent values of a series of correlated observations.
; Specify a critical significance level to test the lag-one auto-correlation coefficient.
; monthly climatolgy perharps not enough data points to calculate
;=================================================================================
; sigr = 0.01 or 0.05 If prob < sigr, then the null hypothesis (means are from the same population)
; is rejected and the alternative hypothesis is accepted.
;================================================================================
sigr = 0.05 ;all areas whose value > 95
sX = new(dimsizes(aveX),typeof(aveX),aveX at _FillValue)
printVarSummary(sX)
;hist_anom1 = dim_avg_n_Wrap(hist_anom,1)
;printVarSummary(hist_anom1) ; (time, lat,lon)
sX = equiv_sample_size (hist_anom(lat|:,lon|:,ens|:0,time|:), sigr,0) ; (time,lat,lon)
copy_VarMeta(aveX,sX)
printMinMax(sX,0)
printVarSummary(sX) ; (lat,lon)
;print(sX)
sY = new(dimsizes(aveX),typeof(aveX),aveX at _FillValue)
printVarSummary(sY)
;rcp45_anom1 = dim_avg_n_Wrap(rcp45_anom,1)
;printVarSummary(rcp45_anom1) ; (time, lat,lon)
sY = equiv_sample_size (rcp45_anom(lat|:,lon|:,ens|:0,time|:), sigr,0) ; (time,lat,lon)
copy_VarMeta(aveY,sY)
printMinMax(sY,0)
printVarSummary(sY) ; (lat,lon)
;print(sY)
;sX = dim_num(.not.ismissing(sX )) ;note that this step looks for the number of actual data points,and doesn't assume that all elements contain data
;sY = dim_num(.not.ismissing(sY ) )
;sX=dim_num_n(.not.ismissing(sX(:,:)),0) ; count missing values at each location
;==================================================================================
; Use "ttest" to compute the probabilities.The returned probability is two-tailed.
;==================================================================================
iflag= True ; population variance similar?
prob = new(dimsizes(varY),typeof(varY),varY at _FillValue)
printVarSummary(prob)
;==================================================================================
; A significance of 0.05 returned by ttest would yield 95% for alpha.
;==================================================================================
;prob = ttest(aveX,varX,sX, aveY,varY,sY, iflag, False)
prob = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, iflag, False))
copy_VarCoords(aveY, prob)
prob at long_name = "Monthly Probablities"
printVarSummary(prob) ; (lat,lon)
printMinMax(prob,0) ; min=0 max=100
;print(prob)
;==================================================================================
;Use "Ftest" to compute the significance of variances
;==================================================================================
f_prob = new(dimsizes(varY),typeof(varY),varY at _FillValue)
f_prob = 100.*(1. - ftest(varX,sX,varY,sY,0))
;f_prob = ftest(varX,sX,varY,sY,0)
copy_VarCoords(aveY, f_prob)
printVarSummary(f_prob) ; (lat,lon)
printMinMax(f_prob,0) ; min=0 max=100 min=3.55992e-05 max=1.2968
;print(f_prob)
print("============")
;===================================================================
; Compute difference Hist vs RCP45
;===================================================================
mean_diff=aveX
mean_diff=aveY-aveX
printVarSummary(mean_diff)
var_diff=varX
var_diff=varY-varX
printVarSummary(var_diff)
skew_diff=xSkwMonEns
skew_diff=xSkwMonEns2-xSkwMonEns
printVarSummary(skew_diff)
kurt_diff=xKurMonEns
kurt_diff=xKurMonEns2-xKurMonEns
printVarSummary(kurt_diff)
print("============")
;==================================================================
; Here is a section of my code that draws a color fill plot
;====================================================================
wks = gsn_open_wks("pdf","Hist_RCP45_difference_statistics")
plot = new(4,graphic)
res = True
res at cnFillOn = True ; turn on color
res at gsnDraw = False ; don't draw
res at gsnFrame = False ; don't advance frame
res at cnInfoLabelOn = False ; turn off cn info label
res at cnFillOn = True ; turn on color
res at mpDataBaseVersion = "MediumRes"
res at cnLinesOn = False
res at tiMainString = " "
;res at tiMainOffsetYF = 0.02 ; Move title up a little
res at gsnCenterString = " "
res at gsnLeftString = " "
res at gsnRightString = " "
;res at gsnSpreadColors = True ; use full range of colormap
;res at gsnSpreadColorStart = 2 ; start at color 2
;res at gsnSpreadColorEnd = -3 ; don't use added gray
res at mpOutlineOn = True
res at mpPerimOn = False
res at mpProjection = "Robinson"
res at lbLabelBarOn = True ; turn off individual cb's
res at lbLabelFontHeightF = 0.013
res at lbBoxEndCapStyle = "TriangleBothEnds" ; Added in NCL V6.4.0
res at lbLabelAutoStride = True ; Control labelbar spacing
;res at tiMainFont = "Helvetica" ; Font for title
;res at tiXAxisFont = "Helvetica" ; Font for X axis label
;res at tiYAxisFont = "Helvetica" ; Font for Y axis label
res at gsnMaximize = True ; large format in landscape
;====================================================================
; to have a common label bar, both plots should be set to the same interval
; b/c the label bar is drawn by default from the interval of the first plot.
; Plot1: Mean difference, and overlay ttest
;==========================================================================
res at cnLevelSelectionMode = "ExplicitLevels"
;res at cnLevels = (/-0.5,-0.3,-0.2,-0.1,0,0.1,0.3,0.5/)
;res at cnFillPalette = "precip_diff_12lev"
res at tiMainString = "a) Ensemble mean difference (RCP45 - Hist)"
;res at lbTitleString ="cm" ;colorbar title
;res at lbTitlePosition ="Right" ; title position
res at lbTitleFontHeightF = .013 ; make title smaller
plot(0) = gsn_csm_contour_map(wks,mean_diff,res)
;res at cnFillPalette = "MPL_Blues"
delete(res at cnLevels)
;res at cnLevels = (/0.06,0.07,0.08,0.09,0.10,0.12,0.15,0.17,0.18,0.19,0.2/)
res at tiMainString = "b) Ensemble Variance difference (RCP45 - Hist)"
;res at lbTitleString ="cm" ;colorbar title
;res at lbTitlePosition ="Right" ; title position
res at lbTitleFontHeightF = .013 ; make title smaller
plot(1) = gsn_csm_contour_map(wks,var_diff,res)
;=====================================================================
; Second plot resources and overlay onto first 2 subplots for mean and variance
;=====================================================================
res2 = True ; res2 probability plots
res2 at gsnDraw = False ; Do not draw plot
res2 at gsnFrame = False ; Do not advance frome
res2 at cnLevelSelectionMode = "ExplicitLevels" ; set explicit cnlev
res2 at cnLevels = (/10./) ; only have 1 contour level
res2 at cnInfoLabelOn = False
res2 at cnLinesOn = False ; do not draw contour lines
res2 at cnLineLabelsOn = False ; do not draw contour labels
res2 at cnFillScaleF = 0.6 ; add extra density
res2 at gsnRightString = "" ; Turn off subtitles
res2 at gsnLeftString = ""
res2 at gsnAddCyclic = False ; add cyclic point
res2 at gsnStringFontHeightF = 0.025
;==================================Create arrays to hold series of ttests or make a loop
neof=2
plot2 = new(neof,graphic) ; create graphic array; only needed if paneling
;res2 at gsnCenterString = "a) Ensemble mean difference (RCP45 - Hist)" ; makes nicer subtitles
plot2(0)=gsn_csm_contour(wks,prob,res2)
;res2 at gsnCenterString = "b) Ensemble Variance difference (RCP45 - Hist)"
plot2(1)=gsn_csm_contour(wks,f_prob,res2)
printVarSummary(plot2)
;========================================================================
; 2b: SHADING
; use pattern fill #2 to fill all areas less than the first contour
; less than or equal to 0.1, and use pattern fill #17 to fill all areas greater
; than the first contour greater than or equal to 1.
;
opt = True ; set up parameters for pattern fill
opt at gsnShadeFillType = "pattern" ; specify pattern fill
opt at gsnShadeHigh = 17 ; stipple pattern
opt at gsnShadeDotSizeF = 0.002 ; make dots larger
;==========================================================================
do mm=0, neof -1
plot2(mm) = gsn_contour_shade(plot2(mm),1,10., opt) ; stipple all areas >= 95% contour ; with sigl=0.05:: no gridbox significant!!!
;plot2(mm) = gsn_contour_shade(plot2(mm),0.1,1., opt)
printVarSummary(plot2)
overlay (plot (mm), plot2 (mm))
end do
;==========================================================================
; last 2 subplots (Skewness and Kurtosis difference)
;==========================================================================
delete(res at cnLevels)
;res at cnLevels = (/-2,-1.8,-1.6,-1.4,-1,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1,1.4,1.6,1.8,2/)
res at cnFillPalette = "BlueDarkRed18"
res at tiMainString = "c) Ensemble Shewness difference (RCP45 - Hist)"
;res at lbTitleString ="cm" ;colorbar title
;res at lbTitlePosition ="Right" ; title position
res at lbTitleFontHeightF = .013 ; make title smaller
plot(2) = gsn_csm_contour_map(wks,skew_diff,res)
delete(res at cnLevels)
;res at cnLevels = (/-2,-1.8,-1.6,-1.4,-1,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1,1.4,1.6,1.8,2/)
res at cnFillPalette = "hotcolr_19lev"
res at tiMainString = "c) Ensemble Kurotsis difference (RCP45 - Hist)"
;res at lbTitleString ="cm" ;colorbar title
;res at lbTitlePosition ="Right" ; title position
res at lbTitleFontHeightF = .013 ; make title smaller
plot(3) = gsn_csm_contour_map(wks,kurt_diff,res)
printVarSummary(plot)
;==================================
; panel resources for the 4 plots
;==================================
pres1 = True ; modify the panel plo
pres1 at gsnMaximize = True ; large format in landscape
gsn_panel(wks,plot,(/2,2/),pres1)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: anomaly10-20_PDF.pdf
Type: application/pdf
Size: 78348 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0011.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: anomaly_all_30N-90N_PDF.pdf
Type: application/pdf
Size: 83558 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0012.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: anomaly_all_30S-30N_PDF.pdf
Type: application/pdf
Size: 83014 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/c6ec23ee/attachment-0013.pdf>
More information about the ncl-talk
mailing list