[ncl-talk] Fwd: ttest: probabilities with very small values

Sri.durgesh Nandini-Weiss sri.durgesh.nandini-weiss at uni-hamburg.de
Thu Nov 14 06:30:48 MST 2019


Hello everyone,

Sorry for re-posting the question:

Im am using NCL to analyses and plot changes in maximum sea level height 
anomalies for different periods.

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...)

I understand that if the sample means in two groups are identical, the 
p-values of a t-test is 1.

So in this case should i do the ttest on monthly anomalies or perhaps 
just de-trending the data is sufficient?

The idea is to first remove the seasonal cycle to otherwise probably not 
Gaussian, and than perform the ttest.

Thanx

Sri



-------- Forwarded Message --------
Subject: 	[ncl-talk] ttest: probabilities with very small values
Date: 	Wed, 13 Nov 2019 11:51:45 +0100
From: 	Sri.durgesh Nandini-Weiss via ncl-talk <ncl-talk at ucar.edu>
Reply-To: 	Sri.durgesh Nandini-Weiss 
<sri.durgesh.nandini-weiss at uni-hamburg.de>
To: 	ncl-talk at ucar.edu



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


_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191114/84da5d85/attachment.html>


More information about the ncl-talk mailing list