<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<p>Hello</p>
<p>Performing ttest on the 2 original samples yields the same
result, with no warning nor error.</p>
<p>The original samples mean:</p>
<p>(0) min=2.11077 max=2.11077<br>
(0) min=0.959119 max=0.959119<br>
</p>
<p>and now the p value is (0) min=2.29254e-10 max=2.29254e-10<br>
</p>
<p>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.</p>
<p>I have attached my script and the pdf plot here when using the
(1) original sample+ pdf and when using the (2) anomalies+pdf.</p>
<p>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.<br>
</p>
<p>Would be grateful for any help.<br>
</p>
<p>Sri</p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<p><br>
</p>
<div class="moz-cite-prefix">On 11/14/19 2:43 PM, Dennis Shea wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAOF1d_6b0EHvsTC+JfS7taHaL64d-yFLJwvUBjXTirL4DWoCsw@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<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" moz-do-not-send="true"><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" moz-do-not-send="true"><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"
moz-do-not-send="true"><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"
moz-do-not-send="true">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" moz-do-not-send="true">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" moz-do-not-send="true">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" moz-do-not-send="true">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"
moz-do-not-send="true">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" moz-do-not-send="true">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote>
</div>
</blockquote>
<pre class="moz-signature" cols="72">--
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</pre>
</body>
</html>