# [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
; =============================================================

hist_anom    = f->hist_anom
printVarSummary(hist_anom) ;[time | 240] x [ens | 100] x [lat | 45] x
[lon | 90]

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

```