[ncl-talk] ttest: probabilities with very small values
Dennis Shea
shea at ucar.edu
Thu Nov 14 06:43:36 MST 2019
I have not looked at the script carefully. Some comments
[1] NCL's *ttest*
<https://www.ncl.ucar.edu/Document/Functions/Built-in/ttest.shtml> is
essentially the same as that in *Numerical Recipes*
<https://websites.pmc.ucsc.edu/~fnimmo/eart290c_17/NumericalRecipesinF77.pdf>
[*pg 610; Chapter 14*].
[2] re: warning:*ttest:* encountered 6 cases where s1 and/or s2 were less
than 2.
s1 and s2 refer to sample sizes
[3] *equiv_sample_size*
<https://www.ncl.ucar.edu/Document/Functions/Built-in/equiv_sample_size.shtml>
Apparently, this is returning very small sample sizes* ?*due to very
high [near 1.0] autocorrelation
*?*
For 'fun', try running the script using the full sample size
On Wed, Nov 13, 2019 at 3:51 AM Sri.durgesh Nandini-Weiss via ncl-talk <
ncl-talk at ucar.edu> wrote:
> 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@
> _FillValue)
> varX = new((/nlat,mlon/), typeof(hist_anom), hist_anom@
> _FillValue)
> xSkwMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom@
> _FillValue)
> xKurMonEns = new((/nlat,mlon/), typeof(hist_anom), hist_anom@
> _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/619ad74a/attachment.html>
More information about the ncl-talk
mailing list