[ncl-talk] Sub-sampling

[Dennis Shea]

https://www.ncl.ucar.edu/Document/Functions/Contributed/generate_sample_indices.shtml

-----
General Comments:


A good programming practice is to not 'hard-code' constants.

model = (/.../)
print(model)

M            = dimsizes(model)   ; M=36 ... not 39
nmodels = ispan(1,M,1)          ; (/1,2,3,...,34,35,36/)

samplesize = (/5,10,15,20,25,30,35/)
N            = dimsizes(samplesize)

Smean = new(N,"float")
            :
dPDJF = asciiread("pr_ens_percent_change_boxplot_DJF.txt",M,"float")
             :
do ni=0,N-1

------
???

Should not
     nmodels = ispan(1,M,1)
be
    nmodels = ispan(0,M-1,1)     ; (/0,1,...,34/)
----

Also, after you upgrade to 6.3.0 or 6.3.1,  rather than

   sdPDJF = dPDJF(generate_sample_indices(samplesize(ni),0))
                 :
   delete(sdPDJF)

You can use the := syntax which was introduced in 6.1.1

   http://www.ncl.ucar.edu/prev_releases.shtml#6.1.1

   sdPDJF := dPDJF(generate_sample_indices(samplesize(ni),0))


------

I am not sure what you mean by:

"I just have a question about how I can do 50 iterations of this to
obtain different models for each time I sub-sample."

Your sample size is 36 ... Not sure how you can sub sample this 50
times and get "different models".

---

If you specify a sample size of (say) 5 and you want each sub sample
of size 5 to have unique indices, I would suggest:

  M  = dimsizes(model)       ; # models
  N  = 5   ; sample size
  mu = generate_unique_indices(M)

  NU = M/N   ; 7

   nuStrt = 0
   nuLast= N-1
   do nu=0,NU-1
     im = mu(nuStrt:nuLast)    ; N unique subscripts
          :
   end do

In the bove the last 'mu' [ ie:  mu(M-1) ] would not be included.

Perhaps, I misunderstand?

What do you have against of sampling with replacement?

On Wed, Oct 28, 2015 at 6:41 AM, Melissa Lazenby <M.Lazenby at sussex.ac.uk> wrote:
> Hi
>
> I have managed to sub-sample 39 models using the generate_sample_indices
> function.
>
> I just have a question about how I can do 50 iterations of this to obtain
> different models for each time I sub-sample.
>
> I think I need to use a do loop but I am not sure as to what the syntax
> should be in my code as I have a loop already to loop through the different
> number of models I sub-sample.
>
> Below is my code currently that is working for just 1 iteration.
>
> I have commented out the attempt at adding the addition iterations loop.
>
> Many thanks!
>
> Kind Regards
> Melissa
>
> ; ==============================================================
> ; sub-sampling_1.ncl
> ; ==============================================================
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
> load "./generate_sample_indices.ncl"
>
> begin
>
> model = (/"ACCESS1-0", "ACCESS1-3", "bcc-csm1-1", "CanESM2", "CCSM4",
> "CESM1-BGC", "CESM1-CAM5", "CMCC-CM", "CMCC-CMS", "CNRM-CM5",
> "CSIRO-Mk3-6-0", "EC-EARTH", "FGOALS-g2", "FIO-ESM","GFDL-CM3",
> "GFDL-ESM2G", "GFDL-ESM2M", "GISS-E2-H_p1", "GISS-E2-H_p2", "GISS-E2-H_p3",
> "GISS-E2-R_p1", "GISS-E2-R_p2", "GISS-E2-R_p3", "HadGEM2-AO", "HadGEM2-CC",
> "HadGEM2-ES", "inmcm4", "IPSL-CM5A-LR", "IPSL-CM5A-MR",
> "IPSL-CM5B-LR","MIROC5", "MPI-ESM-LR", "MPI-ESM-MR", "MRI-CGCM3",
> "NorESM1-ME", "NorESM1-M"/)
>
> nmodels =
> (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39/)
>
> samplesize = (/5,10,15,20,25,30,35/)
> ntimes = 50
>
>
> printVarSummary(samplesize)
>
> Smean = new(7,"float")
> Stdand = new(7,"float")
> Srange = new(7,"float")
> S25 = new(7,"float")
> S75 = new(7,"float")
> Smedian = new(7,"float")
> ;stat = new(7,"float")
>
> dPDJF = asciiread("pr_ens_percent_change_boxplot_DJF.txt",39,"float")
>
>  printVarSummary(dPDJF)
>  print(dPDJF)     ; Print the values
>
>         ;do it=0, ntimes-1
>          do ni=0,6
>
>             sdPDJF = dPDJF(generate_sample_indices(samplesize(ni),0))
>
>               dimt = dimsizes(sdPDJF)   ; should be 5,10,15,20,25,30,35
>                x25  = round(.25*dimt,3)-1     ; -1 to account for NCL
> indexing starting
>                x75  = round(.75*dimt,3)-1     ; at 0
>
>          printVarSummary(sdPDJF)
>
>       ;opt = True
>       ;opt at PrintStat = True
>
>          ;stat(ni) = stat_dispersion(sdPDJF,opt)
>
>          Smean(ni) = avg(sdPDJF)
>          Stdand(ni) = stddev(sdPDJF)
>          Srange(ni) = max(sdPDJF)-min(sdPDJF)
>          S25(ni) = sdPDJF(x25)
>          S75(ni) = sdPDJF(x75)
>          Smedian(ni) = dim_median(sdPDJF)
>
>
>         print(Smean)
>         print(Stdand)
>         print(Srange)
>         print(S25)
>         print(S75)
>         print(Smedian)
>
>          delete(sdPDJF)
>
>          end do
>         ;end do
> end
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>