<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<p><br>
</p>
<div class="moz-forward-container">Hello everyone,</div>
<div class="moz-forward-container"><br>
</div>
<div class="moz-forward-container">Sorry for re-posting the
question:</div>
<div class="moz-forward-container"><br>
</div>
<div class="moz-forward-container">
<p>Im am using NCL to analyses and plot changes in maximum sea
level height anomalies for different periods.</p>
<p>However, when i compute t-test, the p-value is equal to 1
because my computed mean values (on the anomalies) are very very
small, almost 0 (which is what i should expect...)<br>
</p>
<p>I understand that if the sample means in two groups are
identical, the p-values of a t-test is 1.</p>
<p>So in this case should i do the ttest on monthly anomalies or
perhaps just de-trending the data is sufficient?</p>
<p>The idea is to first remove the seasonal cycle to otherwise
probably not Gaussian, and than perform the ttest.</p>
<p>Thanx</p>
<p>Sri<br>
</p>
</div>
<div class="moz-forward-container"><br>
</div>
<div class="moz-forward-container"><br>
</div>
<div class="moz-forward-container">-------- Forwarded Message
--------
<table class="moz-email-headers-table" cellspacing="0"
cellpadding="0" border="0">
<tbody>
<tr>
<th valign="BASELINE" nowrap="nowrap" align="RIGHT">Subject:
</th>
<td>[ncl-talk] ttest: probabilities with very small values</td>
</tr>
<tr>
<th valign="BASELINE" nowrap="nowrap" align="RIGHT">Date: </th>
<td>Wed, 13 Nov 2019 11:51:45 +0100</td>
</tr>
<tr>
<th valign="BASELINE" nowrap="nowrap" align="RIGHT">From: </th>
<td>Sri.durgesh Nandini-Weiss via ncl-talk
<a class="moz-txt-link-rfc2396E" href="mailto:ncl-talk@ucar.edu"><ncl-talk@ucar.edu></a></td>
</tr>
<tr>
<th valign="BASELINE" nowrap="nowrap" align="RIGHT">Reply-To:
</th>
<td>Sri.durgesh Nandini-Weiss
<a class="moz-txt-link-rfc2396E" href="mailto:sri.durgesh.nandini-weiss@uni-hamburg.de"><sri.durgesh.nandini-weiss@uni-hamburg.de></a></td>
</tr>
<tr>
<th valign="BASELINE" nowrap="nowrap" align="RIGHT">To: </th>
<td><a class="moz-txt-link-abbreviated" href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a></td>
</tr>
</tbody>
</table>
<br>
<br>
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
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 only one warning:<br>
<br>
warning:ttest: encountered 6 cases where s1 and/or s2 were less
than 2. Output set to missing values in these cases.<br>
<br>
Heres my script:<br>
<br>
;========================================================================================<br>
; Concepts illustrated in my script:
<a class="moz-txt-link-abbreviated" href="mailto:sri.durgesh.nandini-weiss@uni-hamburg.de">sri.durgesh.nandini-weiss@uni-hamburg.de</a><br>
; The attached calculates the ensemble (N=100) mean wave height,
standard deviation, skewness and kurtosis between 1981-2006 and
2080-2100<br>
; Then it plots ensembl mean higher moments in statistics by
computing 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 ("anom_hist_zo.nc", "r")<br>
hist_anom = f->hist_anom<br>
printVarSummary(hist_anom) ;[time | 240] x [ens | 100] x [lat |
45] x [lon | 90]<br>
<br>
f1 = addfile ("anom_rcp45_zo.nc", "r")<br>
rcp45_anom = f1->rcp45_anom<br>
printVarSummary(rcp45_anom) ;[time | 240] x [ens | 100] x [lat
| 45] 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 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 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),
rcp45_anom@_FillValue)<br>
varY = new((/nlat,mlon/), typeof(rcp45_anom),
rcp45_anom@_FillValue)<br>
xSkwMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom),
rcp45_anom@_FillValue)<br>
xKurMonEns2 = new((/nlat,mlon/), typeof(rcp45_anom),
rcp45_anom@_FillValue)<br>
<br>
do nmo=0,nmos-1<br>
work2 := reshape(rcp45_anom(nmo::nmos,:,:,:)
,(/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 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 series of correlated observations.<br>
; Specify a critical significance level to test the lag-one
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 are from the same population)<br>
; is rejected and the alternative hypothesis is accepted.<br>
;================================================================================<br>
<br>
sigr = 0.05 ;all areas
whose 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|:),
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|:), 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 two-tailed.<br>
;==================================================================================<br>
<br>
iflag= True ;
population 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 class="moz-txt-link-abbreviated" href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a class="moz-txt-link-freetext" href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></div>
</body>
</html>