<div dir="ltr"><div>I have not looked at the script carefully. Some comments</div><div><br></div><div>[1] NCL's <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml" target="_blank"><b>ttest</b></a> is essentially the same as that in <a href="https://websites.pmc.ucsc.edu/~fnimmo/eart290c_17/NumericalRecipesinF77.pdf" target="_blank"><b>Numerical Recipes</b></a> [<b>pg 610; Chapter 14</b>].</div><div><div>[2] re: warning:<b>ttest:</b> encountered 6 cases where s1 and/or s2 were less than 2. <br></div><div> s1 and s2 refer to sample sizes<b> <br></b></div></div><div>[3] <a href="https://www.ncl.ucar.edu/Document/Functions/Built-in/equiv_sample_size.shtml"><b>equiv_sample_size</b></a></div><div><b> </b>Apparently, this is returning very small sample sizes<b> ?</b>due to very high [near 1.0] autocorrelation<b>?<br></b></div><div><b> </b>For 'fun', try running the script<b> </b>using the full sample size<b><br></b></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Nov 13, 2019 at 3:51 AM Sri.durgesh Nandini-Weiss via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Good morning everyone,<br>
<br>
I was wondering to get some advice on the below.<br>
<br>
I have a working script for computing ttest, which is computed on <br>
anomalies, however, the values are incredibly small-<br>
<br>
at the significance of 0.05 returned by ttest would yield 95% for alpha.<br>
<br>
Probablities : min=0 max=0.000435114<br>
SO im wondering if it actually makes sense to try to plot it. I have <br>
only one warning:<br>
<br>
warning:ttest: encountered 6 cases where s1 and/or s2 were less than 2. <br>
Output set to missing values in these cases.<br>
<br>
Heres my script:<br>
<br>
;========================================================================================<br>
; Concepts illustrated in my script: <br>
<a href="mailto:sri.durgesh.nandini-weiss@uni-hamburg.de" target="_blank">sri.durgesh.nandini-weiss@uni-hamburg.de</a><br>
; The attached calculates the ensemble (N=100) mean wave height, <br>
standard deviation, skewness and kurtosis between 1981-2006 and 2080-2100<br>
; Then it plots ensembl mean higher moments in statistics by computing <br>
20yrs anomalies 4 (3) moments in panel:<br>
; (a) Monthly climatology and monthly anomalies<br>
; (b) Ensemble spread (std. deviation)<br>
; (c) skewness and (d) kurtosis<br>
<br>
;==============================================================<br>
; Open the file: Read only the user specified period<br>
; =============================================================<br>
<br>
f = addfile ("<a href="http://anom_hist_zo.nc" rel="noreferrer" target="_blank">anom_hist_zo.nc</a>", "r")<br>
hist_anom = f->hist_anom<br>
printVarSummary(hist_anom) ;[time | 240] x [ens | 100] x [lat | 45] x <br>
[lon | 90]<br>
<br>
f1 = addfile ("<a href="http://anom_rcp45_zo.nc" rel="noreferrer" target="_blank">anom_rcp45_zo.nc</a>", "r")<br>
rcp45_anom = f1->rcp45_anom<br>
printVarSummary(rcp45_anom) ;[time | 240] x [ens | 100] x [lat | 45] <br>
x [lon | 90]<br>
<br>
<br>
dimx = dimsizes(hist_anom)<br>
ntim = dimx(0) ; 240<br>
nens = dimx(1) ; 100<br>
nlat = dimx(2) ; 45<br>
mlon = dimx(3) ; 90<br>
<br>
nmos = 12<br>
nyrs = ntim/nmos ; 20<br>
printVarSummary(nyrs)<br>
<br>
;==================================================================<br>
; Compute first four moments (average, variance, skewness, and kurtosis)<br>
; over the time and ensemble dimension together: for Hist, and then RCP45<br>
;==================================================================<br>
<br>
aveX = new((/nlat,mlon/), typeof(hist_anom), hist_anom@_FillValue)<br>
varX = new((/nlat,mlon/), typeof(hist_anom), hist_anom@_FillValue)<br>
xSkwMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom@_FillValue)<br>
xKurMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom@_FillValue)<br>
<br>
do nmo=0,nmos-1<br>
work := reshape(hist_anom(nmo::nmos,:,:,:) ,(/nyrs*nens,nlat,mlon/))<br>
if (nmo.eq.0) then<br>
printVarSummary(work) ; (2000,45,90)<br>
end if<br>
xStat = dim_stat4_n(work,0 )<br>
end do<br>
printVarSummary(xStat) ; (4,nlat,mlon)<br>
aveX = xStat(0,:,:) ; (nlat, mlon)<br>
varX = xStat(1,:,:)<br>
xSkwMonEns = xStat(2,:,:)<br>
xKurMonEns = xStat(3,:,:)<br>
<br>
copy_VarCoords(hist_anom(0,0,:,:), aveX)<br>
aveX@long_name = "Monthly Climatology (Average)"<br>
aveX@LONG_NAME = "Monthly Climatology over all Ensemble Members"<br>
<br>
copy_VarCoords(hist_anom(0,0,:,:), varX)<br>
varX@long_name = "Monthly Interannual Variability (sample variance)"<br>
varX@LONG_NAME = "Monthly Interannual Variability over all Ensemble <br>
Members"<br>
<br>
copy_VarCoords(hist_anom(0,0,:,:), xSkwMonEns)<br>
xSkwMonEns@long_name = "Monthly Skewness"<br>
copy_VarCoords(hist_anom(0,0,:,:), xKurMonEns)<br>
xKurMonEns@long_name = "Monthly Kurtosis"<br>
<br>
print("============")<br>
<br>
;================================================================================<br>
; dim average on month dimension; does this change the plotted results <br>
and prob?<br>
;================================================================================<br>
<br>
printVarSummary(aveX) ; (lat,lon)<br>
printMinMax(aveX,0)<br>
printVarSummary(varX) ; (lat,lon)<br>
printMinMax(varX,0)<br>
printVarSummary(xSkwMonEns) ; (lat,lon)<br>
printMinMax(xSkwMonEns,0)<br>
printVarSummary(xKurMonEns) ; (lat,lon)<br>
printMinMax(xKurMonEns,0)<br>
print("============")<br>
<br>
;==================================================================<br>
; Now repeat for RCP45<br>
;==================================================================<br>
<br>
<br>
aveY = new((/nlat,mlon/), typeof(rcp45_anom), <br>
rcp45_anom@_FillValue)<br>
varY = new((/nlat,mlon/), typeof(rcp45_anom), <br>
rcp45_anom@_FillValue)<br>
xSkwMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom), <br>
rcp45_anom@_FillValue)<br>
xKurMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom), <br>
rcp45_anom@_FillValue)<br>
<br>
do nmo=0,nmos-1<br>
work2 := reshape(rcp45_anom(nmo::nmos,:,:,:) <br>
,(/nyrs*nens,nlat,mlon/))<br>
if (nmo.eq.0) then<br>
printVarSummary(work2) ; (2000,45,90)<br>
end if<br>
xStat2 = dim_stat4_n(work2,0 )<br>
end do<br>
printVarSummary(xStat2) ; (4,nlat,mlon)<br>
aveY = xStat2(0,:,:) ; (nlat, mlon)<br>
varY = xStat2(1,:,:)<br>
xSkwMonEns2 = xStat2(2,:,:)<br>
xKurMonEns2 = xStat2(3,:,:)<br>
<br>
copy_VarCoords(rcp45_anom(0,0,:,:), aveY)<br>
aveY@long_name = "Monthly Climatology (Average)"<br>
aveY@LONG_NAME = "Monthly Climatology over all Ensemble Members"<br>
<br>
copy_VarCoords(rcp45_anom(0,0,:,:), varY)<br>
varY@long_name = "Monthly Interannual Variability (sample variance)"<br>
varY@LONG_NAME = "Monthly Interannual Variability over all Ensemble <br>
Members"<br>
<br>
copy_VarCoords(rcp45_anom(0,0,:,:), xSkwMonEns2)<br>
xSkwMonEns2@long_name = "Monthly Skewness"<br>
copy_VarCoords(rcp45_anom(0,0,:,:), xKurMonEns2)<br>
xKurMonEns2@long_name = "Monthly Kurtosis"<br>
<br>
print("============")<br>
<br>
printVarSummary(aveY) ; (lat,lon)<br>
printMinMax(aveY,0)<br>
printVarSummary(varY) ; (lat,lon)<br>
printMinMax(varY,0)<br>
printVarSummary(xSkwMonEns2) ; (lat,lon)<br>
printMinMax(xSkwMonEns2,0)<br>
printVarSummary(xKurMonEns2) ; (lat,lon)<br>
printMinMax(xKurMonEns2,0)<br>
print("============")<br>
<br>
;================================================================================<br>
; equiv_sample_size: Estimates the number of independent values of a <br>
series of correlated observations.<br>
; Specify a critical significance level to test the lag-one <br>
auto-correlation coefficient.<br>
; monthly climatolgy perharps not enough data points to calculate<br>
;=================================================================================<br>
; sigr = 0.01 or 0.05 If prob < sigr, then the null hypothesis (means <br>
are from the same population)<br>
; is rejected and the alternative hypothesis is accepted.<br>
;================================================================================<br>
<br>
sigr = 0.05 ;all areas whose <br>
value > 95<br>
<br>
sX = new(dimsizes(aveX),typeof(aveX),aveX@_FillValue)<br>
printVarSummary(sX)<br>
<br>
sX = equiv_sample_size (hist_anom(lat|:,lon|:,ens|:0,time|:), <br>
sigr,0) ; (time,lat,lon)<br>
copy_VarMeta(aveX,sX)<br>
printMinMax(sX,0)<br>
printVarSummary(sX) ; (lat,lon)<br>
;print(sX)<br>
<br>
sY = new(dimsizes(aveX),typeof(aveX),aveX@_FillValue)<br>
printVarSummary(sY)<br>
<br>
sY = equiv_sample_size (rcp45_anom(lat|:,lon|:,ens|:0,time|:), <br>
sigr,0) ; (time,lat,lon)<br>
copy_VarMeta(aveY,sY)<br>
printMinMax(sY,0)<br>
printVarSummary(sY) ; (lat,lon)<br>
;print(sY)<br>
<br>
<br>
;==================================================================================<br>
; Use "ttest" to compute the probabilities.The returned probability is <br>
two-tailed.<br>
;==================================================================================<br>
<br>
iflag= True ; population <br>
variance similar?<br>
prob = new(dimsizes(varY),typeof(varY),varY@_FillValue)<br>
printVarSummary(prob)<br>
<br>
;==================================================================================<br>
; A significance of 0.05 returned by ttest would yield 95% for alpha.<br>
;==================================================================================<br>
<br>
prob = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, iflag, False))<br>
copy_VarCoords(aveY, prob)<br>
prob@long_name = " Probablities"<br>
printVarSummary(prob) ; (lat,lon)<br>
printMinMax(prob,0) ; min=0 max=100<br>
;print(prob)<br>
<br>
;==================================================================================<br>
<br>
Best<br>
<br>
Sri<br>
<br>
<br>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>