[ncl-talk] ttest: probabilities with very small values
Sri.durgesh Nandini-Weiss
sri.durgesh.nandini-weiss at uni-hamburg.de
Wed Nov 13 03:51:45 MST 2019
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
; 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) ;[time | 240] x [ens | 100] x [lat | 45] x
[lon | 90]
f1 = addfile ("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
More information about the ncl-talk
mailing list