scilab/modules/statistics/macros/histc.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2017 - 2020 - Samuel GOUGEON
   3 //
   4 // This file is hereby licensed under the terms of the GNU GPL v2.0,
   5 // pursuant to article 5.3.4 of the CeCILL v.2.1.
   6 // This file was originally licensed under the terms of the CeCILL v2.1,
   7 // and continues to be available under such terms.
   8 // For more information, see the COPYING file which you should have received
   9 // along with this program.
  10
  11
  12 function [heights, jokers, binsOut, ind] = histc(data, bins, options)
  13 // Scilab 5.5:
  14 //  1) histc(bins, data, normalization)
  15 //  2) By défaut, histc() worked in "densityNorm,normWith:all" mode
  16 // Scilab 6.0:
  17 //  1) histc(data, bins, options)
  18 //  2) By default, histc() works in "counts" mode
  19
  20 //    histc(data)
  21 //    histc(data, nbins)
  22 //    histc(data, binsWidth<0)
  23 //    histc(data, binsEdges)
  24 //    histc(data, binsValues [, "discrete"])
  25 //    histc(data, binsAlgo)
  26 //    histc(data, ..., options)
  27 //    histc(nbins, data)                        // deprecated
  28 //    histc(binsEdges, data)                    // deprecated
  29 //    histc(binsEdges, data, normalization=%t)  // deprecated
  30 //     heights                        = histc(data ..)
  31 //    [heights, jokers]               = histc(data ..)
  32 //    [heights, jokers, binsOut]      = histc(data ..)
  33 //    [heights, jokers, binsOut, ind] = histc(data ..)
  34 //
  35 //    Binning algorithm:
  36 //      "sqrt" | "sturges" | "freediac"
  37 //    Binning options:
  38 //      "discrete"  // , "smartbins"
  39 //    Scale options:
  40 //      "counts" (default) | "countsNorm" | "density" | "densityNorm"
  41 //    Normalization "normWith: " options:  Priorities: "all" > "out" > ..
  42 //     decimal, complex numbers: "leftout rightout out inf nan zeros all"
  43 //                     integers: "leftout rightout (out=) all"
  44 //                         Text: "empty out all"
  45 //                  polynomials: "(out=) all"
  46 //       vector or single string with space-separated flags (case-unsensitive)
  47 //
  48 //    Complex numbers are sorted by real parts, then by imaginary parts
  49 //    jokers =
  50 //      * integers: [leftout, rightout] (continuous) or [out] (discrete)
  51 //      * decimal or complex numbers, full or sparse:
  52 //            continuous: [leftout, rightout, zeros, nan, inf]
  53 //              discrete: [out, 0, zeros, nan, inf]
  54 //      * polynomials: [out, 0, 0, nan, inf]
  55 //      * texts: continuous: [leftout, rightout, ""]; discrete: [out, 0, ""]
  56 //
  57
  58     [lhs, rhs] = argn();
  59     fname = "histc"
  60     [heights, jokers, binsOut, ind] = ([],[],[],[])
  61     [Nout, Nleftout, Nzeros, Nrightout, Ninf, Nnan] = (0,0,0,0,0,0)
  62
  63     // PARSING INPUT ARGUMENTS
  64     // =======================
  65     if rhs<1 | rhs>3 then
  66         msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
  67         error(msprintf(msg, fname, 1, 3))
  68     end
  69     // Old 5.5.x mode management
  70     // -------------------------
  71     mode55 = %f
  72     // Check for Scilab 5.5 back-compatibility
  73     if (isdef("data","l") & isdef("bins", "l") & type(bins)>0) & ..
  74         size(data,"*")< size(bins,"*") & or(type(data)==[1 8]) & type(bins)~=10
  75         msg = _("%s: data are now expected in argument #1 => arguments #1 and #2 switched\n")
  76         warning(msprintf(msg, fname))
  77         mode55 = %t
  78         [data, bins] = (bins, data)
  79         tmp = isdef("options","l") && type(options)==4
  80         normarg = %f
  81         norma = %f
  82         if isdef("normalization","l")
  83             normarg = %t
  84             if type(normalization)==4 then
  85                 norma = normalization(1)
  86             end
  87         else
  88             if rhs<3 then
  89                 norma = %t
  90             else
  91                 if isdef("options","l") & type(options)==4
  92                     normarg = %t
  93                     norma = options(1)
  94                 end
  95             end
  96         end
  97         if norma then
  98             options = ["densityNorm" "normWith:all"]
  99         else
 100             options = "counts"
 101         end
 102         if normarg
 103             msg = _("%s: ""normalization"" boolean argument #3 is obsolete. ""%s"" option used.\n")
 104             warning(msprintf(msg, fname, options(1)))
 105         end
 106     end
 107     // -----------------------
 108     if ~isdef("data","l") | type(data)==0 then
 109         msg = _("%s: Argument #%d: %s\n")
 110         error(msprintf(msg, fname, 1, "data array missing"))
 111     end
 112     if ~isdef("bins","l") | type(bins)==0 then
 113         bins = []
 114     end
 115     if ~isdef("options","l") then
 116         options = "counts"
 117     end
 118     // Now, bins, data and options are defined. Let's check their types and sizes
 119
 120     // CHECKING INPUT ARGUMENTS AND SETTING DEFAULTS
 121     // =============================================
 122     // Options
 123     // -------
 124     if options~=[]
 125         if type(options)~=10 then
 126             msg = _("%s: Argument #%d: Text expected.\n")
 127             error(msprintf(msg, fname, 3))
 128         end
 129         options = tokens(strcat(convstr(options),","),",")
 130     end
 131     discrete = grep(options, "discrete")~=[] | type(data)==2 // polynomials
 132
 133     // Data
 134     // ----
 135     if isempty(data) then
 136         return
 137     end
 138     if ~or(type(data)==[1 2 5 8 10]) then
 139         msg = _("%s: Data in argument #%d: Numbers or polynomials or texts expected.\n")
 140         error(msprintf(msg, fname, 1))
 141     end
 142     data0 = data    // Saving
 143     // If encoded integers (including u-int64): => decimal (for dsearch(), etc)
 144     if type(data) == 8 then
 145         data = double(data)
 146     end
 147     // If sparse numbers: we remove zeros
 148     if type(data)==5
 149         [IJ, data] = spget(data)
 150     end
 151     // If there are only %inf or %nan => return ==!!!!!!!!!!!!!!!!!!!!!!!!!!!
 152     if ~discrete & or(type(data0)==[1 5]) then
 153         L = data
 154         tmp = isnan(data)
 155         Nnan = sum(tmp)
 156         L(tmp) = []
 157         tmp = isinf(data)
 158         Ninf = sum(tmp)
 159         L(tmp) = []
 160         if data==[]
 161             msg = _("%s: data have only %inf or/and %nan values.\n")
 162             warning(msprintf(msg, fname))
 163             if type(data0)==1
 164                 Nzeros = sum(data==0)
 165             else
 166                 Nzeros = size(data0,"*") - Nnan - Ninf
 167             end
 168             Nleftout = 0
 169             Nrightout = 0
 170             if length(bins)>1
 171                 b = bins(~isinf(bins))
 172                 if 0<min(b)
 173                     Nleftout = Nzeros
 174                 end
 175                 if 0 > max(b)
 176                     Nrightout = Nzeros
 177                 end
 178             end
 179             jokers = [Ninf Nnan Nzeros Nleftout Nrightout]
 180             return
 181         end
 182     end
 183     // Polynomials: setting Nan and Ninf
 184     if type(data)==2
 185         L = data
 186         tmp = isnan(L)
 187         Nnan = sum(tmp)
 188         L(tmp) = []     // to ensure excluding any intersection
 189         tmp = isinf(L)
 190         Ninf = sum(tmp)
 191     end
 192
 193     // Bins
 194     // ----
 195     // Setting a default binning mode
 196     if bins==[] then
 197         // histc(data): Default bins:
 198         //  * strings: discrete, from unique()
 199         //  * polynomials: discrete, from defBins set hereabove
 200         //  * encoded integers: "sqrt". if strange() < "sqrt" => binsWidth = 1
 201         //  * decimal, complex: binMethod option or default "sqrt"
 202         if type(data)==10       // texts
 203             discrete = %t
 204             bins = unique(data)(:)'
 205             bins = bins(:)'
 206
 207         elseif type(data)==2    // Polynomials
 208             L = data
 209             tmp = isnan(L)
 210             L(tmp) = []
 211             bins = []
 212             while L~=[]
 213                 bins = [bins L(1)]
 214                 L(L==L(1)) = []
 215             end
 216             clear L
 217
 218         else    // decimals, complexes, encoded integers
 219             if ~discrete
 220                 bins = "sqrt"
 221             else
 222                 bins = unique(data)(:).'
 223             end
 224         end
 225     end
 226     //
 227     if ~discrete
 228         if type(bins)==8 | (type(bins)==1 & ~isreal(bins))
 229             bins = real(double(bins))
 230         end
 231         if size(bins,"*")==1 & type(bins)==10 then
 232             // binsAlgo
 233             // --------
 234             if ~or(convstr(bins)==["sqrt" "sturges" "freediac"])
 235                 msg = _("%s: Argument #%d: wrong value for binning algo")
 236                 error(msprintf(msg, fname, 2))
 237             end
 238             // We set the number of bins. Their edges are set later.
 239             Nvalid = size(data,"*") - Nnan - Ninf
 240             select bins
 241             case "sturges"
 242                 bins = ceil(1 + log2(Nvalid))
 243             case "freediac"
 244                 tmp = real(data)
 245                 binWidth = 2*iqr(tmp)* Nvalid^(-1/3)
 246                 if binWidth~=0
 247                     bins = round((max(tmp) - min(tmp)) / binWidth)
 248                 else
 249                     bins = 2
 250                 end
 251             else  // "sqrt"
 252                 bins = max(1, round(sqrt(Nvalid)))
 253             end
 254         end
 255
 256         if length(bins)==1 & type(bins)==1  then
 257             // The number of bins or their width is provided
 258             if type(data)==10   // texts
 259                 msg = _("%s: Argument #2: Please provide bins edges or values or leave choosing default bins.\n")
 260                 error(msprintf(msg, fname))
 261             end
 262             // Case polynomials: not possible here <= discrete previously set to %t
 263             // Cases: Encoded integers, decimal or complex numbers
 264             if bins>=0 & (bins~=floor(bins) | bins==0) | isnan(bins) | bins==%inf
 265                 msg = _("%s: Argument #%d: non-zero decimal integer expected.\n")
 266                 error(msprintf(msg, fname, 2))
 267             elseif bins==-%inf
 268                 msg = _("%s: Argument #%d: decimal number > -Inf expected.\n")
 269                 error(msprintf(msg, fname, 2))
 270             end
 271
 272             // If data are complex-encoded, they will be sorted by increasing
 273             //  real parts, then by increasing imaginary parts.
 274             b = data(~isnan(data) & ~isinf(data))
 275             b = real(b)
 276             mind = min(b)
 277             maxd = max(b)
 278             if type(data0)==1
 279                 if bins > 0
 280                     // Required number of bins
 281                     if (mind == maxd) then
 282                         mind = mind*(1 - 0.5/bins);
 283                         maxd = maxd*(1 + 0.5/bins);
 284                     end
 285                     bins = linspace(mind, maxd, bins+1);          // Bins edges
 286                 else
 287                     // Required bins width
 288                     if (mind == maxd) then
 289                         bins = mind + [1 -1]*bins/2
 290                     else
 291                         // Edges set à 0.0, width = |bins|, from mind- to maxd+
 292                         tmp = -bins
 293                         bins = (floor(mind/tmp):ceil(maxd/tmp))*tmp
 294                     end
 295                 end
 296             else
 297                 // Encoded integers
 298                 if bins > 0
 299                     // Required number of bins
 300                     tmp = maxd - mind + 1
 301                     if tmp <= bins & tmp <= 2^16  // integer-centered bins of width==1
 302                         bins = (mind:(maxd+1)) - 0.5              // Bins edges
 303                     else
 304                         ds = (maxd - mind)/bins/2
 305                         bins = linspace(mind-ds, maxd+ds, bins+1) // Bins edges
 306                     end
 307                 else
 308                     // Required bins width
 309                     tmp = max(round(-bins),1) // Bins edges @ -0.5 modulo |int(bins)|
 310                     bins = (floor(mind/tmp):ceil(maxd/tmp))*tmp - 0.5 // Bins edges
 311                 end
 312             end
 313
 314         else
 315             // Bins edges are provided: we check that edges are sorted
 316             // Texts
 317             // -----
 318             if type(data0)==10
 319                 if type(bins)~=10
 320                     msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
 321                     error(msprintf(msg, fname, 1, 2))
 322                 end
 323                 if ~and(bins==unique(bins))
 324                     msg = _("%s: Argument #%d: Elements must be sorted in increasing order.\n")
 325                     error(msprintf(msg, fname, 1))
 326                 end
 327
 328             // Numbers
 329             // -------
 330             else
 331                 if ~or(type(bins)==[1 8])
 332                     msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
 333                     error(msprintf(msg, fname, 1, 2))
 334                 end
 335                 bins = bins(~isnan(bins))
 336                 if min(diff(bins)) <= 0 then
 337                     msg = _("%s: Argument #%d: Elements must be in increasing order.\n")
 338                     error(msprintf(msg, fname, 1))
 339                 end
 340                 // We manage -%inf open left bin
 341                 if bins(1)==-%inf
 342                     tmp = min(data(~isinf(data)))
 343                     if bins(2)<=tmp
 344                         bins(1) = []
 345                     else
 346                         bins(1) = tmp
 347                     end
 348                 end
 349                 // We manage %inf open right bin
 350                 if bins($)==%inf
 351                     tmp = max(data(~isinf(data)))
 352                     if bins($-1)>=tmp
 353                         bins($) = []
 354                     else
 355                         bins($) = tmp
 356                     end
 357                 end
 358             end
 359         end
 360          // if ~discrete
 361
 362     else    // if discrete
 363         if type(data0)==10
 364             if type(bins)~=10
 365                 msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
 366                 error(msprintf(msg, fname, 1, 2))
 367             end
 368             [bins, dataOrder] = unique(bins)
 369
 370         elseif or(type(data)==[1 8])
 371             if ~or(type(bins)==[1 8])
 372                 msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
 373                 error(msprintf(msg, fname, 1, 2))
 374             end
 375             bins = real(double(bins))
 376             bins(isnan(bins)) = []
 377             [bins, dataOrder] = unique(bins)
 378         end
 379     end
 380     binsOut = bins
 381
 382     // PROCESSING
 383     // ==========
 384     // jokers =
 385     //  * integers: [leftout, rightout] (continuous) or [out] (discrete)
 386     //  * decimal or complex numbers, full or sparse:
 387     //       continuous: [inf, nan, zeros, leftout, rightout]
 388     //         discrete: [inf, nan, zeros, out]
 389     //  * polynomials: [inf, nan, out]
 390     //  * texts: continuous: [empty, leftout, rightout]; discrete: [empty, out]
 391     if type(data0)~=2 then
 392         hmode = "c"            // continuous
 393         if discrete then
 394             hmode = "d"        // discrete
 395         end
 396         if or(type(data0)==[1 5])
 397             if type(data0)==1
 398                 tmp = data0
 399             else
 400                 tmp = data     // sparse
 401             end
 402             tmp1 = isinf(real(tmp))
 403             tmp2 = isinf(imag(tmp))
 404             tmp3 = isnan(tmp) | (tmp1 & tmp2)
 405             if ~discrete & ~isreal(tmp,0)
 406                 data(tmp1 | tmp2) = %inf
 407                 data(tmp3) = %nan
 408             end
 409             clear tmp tmp1 tmp2 tmp3
 410             data = real(data)
 411         end
 412         [ind, heights] = dsearch(data, bins, hmode)
 413                                      // ind: bin's number of each data component
 414
 415         if discrete
 416             // We restore the original order of binsValues:
 417             [tmp, k]  = gsort(dataOrder, "g","i")
 418             [tmp, k2] = gsort(k,"g","i")
 419             bins = bins(k)
 420             binsOut = bins
 421             heights = heights(k)
 422             tmp2 = ind~=0
 423             ind(tmp2) = k2(ind(tmp2))
 424         end
 425
 426     else // Polynomials
 427         // Counts (and location)
 428          ind = zeros(data)
 429          heights = []
 430          for i = 1:size(bins,"*")
 431              tmp = data==bins(i)
 432              heights(1,i) = sum(tmp)
 433              ind(tmp) = i
 434          end
 435          tmp = Ninf - sum(heights(isinf(bins)))     // inf out of the selection
 436          Nout = size(data,"*") - sum(heights) - Nnan - tmp// excludes all nan and inf not in the selection
 437          jokers = [Nout, 0, 0, Nnan, Ninf]
 438     end
 439
 440     // POST-PROCESSING
 441     // ===============
 442     // Counting jokers
 443     // ---------------
 444     if or(type(data0)==[1 5])        // Real and complex numbers, sparse or not
 445         inf = isinf(data)
 446         Ninf = sum(inf)
 447         Nnan = sum(isnan(data))
 448         if type(data0)==1
 449             Nzeros = sum(data==0)
 450         else
 451             Nzeros = size(data0,"*") - nnz(data0)
 452         end
 453         if ~discrete                 // Sparses: 0 is not taken into account
 454             Nleftout  = sum(data(~inf) < bins(1))
 455             Nrightout = sum(data(~inf) > bins($))
 456             jokers = [Nleftout Nrightout Nzeros Nnan Ninf]
 457         else
 458             Nout = size(data0,"*") - sum(heights) - Nnan - Ninf
 459             binsinf = isinf(bins)
 460             if or(binsinf)
 461                 Nout = Nout + sum(heights(binsinf))
 462             end
 463             if type(data0)==5
 464                 Nout = Nout - Nzeros
 465             end
 466             jokers = [Nout 0 Nzeros Nnan Ninf]
 467         end
 468
 469     elseif type(data0)==8            // Encoded integers
 470         if ~discrete
 471             Nleftout  = sum(data < bins(1))
 472             Nrightout = sum(data > bins($))
 473             jokers = [Nleftout Nrightout]
 474             Nout = Nleftout + Nrightout
 475         else
 476             Nout = size(data,"*") - sum(heights)
 477             jokers = Nout
 478         end
 479
 480     elseif type(data0)==10           // Texts
 481         Nempty = sum(data0=="")
 482         Nout = size(data0,"*") - sum(heights)
 483         if ~discrete
 484             tmp = gsort(data,"g","i")
 485             Nleftout = sum(find(strcmp(tmp,bins(1))<0)($)) - Nempty
 486             tmp = sum(find(strcmp(tmp,bins($))>0)(1))
 487             if ~tmp
 488                 Nrightout = 0
 489             else
 490                 Nrightout = size(data,"*") - tmp + 1
 491             end
 492             clear tmp
 493             jokers = [Nleftout, Nrightout, Nempty]
 494         else
 495             if ~or(bins=="")
 496                 Nout = Nout - Nempty
 497             end
 498             jokers = [Nout, 0, Nempty]
 499         end
 500
 501     // Polynomials: jokers already set when setting bins
 502     end
 503
 504     // Tuning the histogram scale
 505     // --------------------------
 506     histmode = "counts"
 507     if discrete then
 508         tmp = grep(options, ["counts" "countsnorm"])
 509     else
 510         tmp = grep(options, ["counts" "countsnorm" "density" "densitynorm"])
 511     end
 512     if tmp~=[] then
 513         histmode = options(tmp($))
 514     end
 515     // jokers =
 516     //  * integers: [leftout, rightout] (continuous) or [out] (discrete)
 517     //  * decimal or complex numbers, full or sparse:
 518     //       continuous: [leftout, rightout, zeros, nan, inf]
 519     //         discrete: [out, 0, zeros, nan, inf]
 520     //  * polynomials: [out, 0, 0, nan, inf].
 521     //  * texts: continuous: [leftout, rightout, empty]; discrete: [out, empty]
 522
 523     // Number of components to take into account:
 524     // Parsing the option "normWith:
 525     //  Continuous: " leftout rightout inf nan zeros"
 526     //    Discrete: " out inf nan zeros"
 527     if grep(histmode, "norm") then
 528         N = sum(heights)
 529         tmp = grep(options, "normwith")
 530         normMode = options(tmp)
 531         if normMode~=[]
 532             gh = "/[\:|\s+]"    // grep() header
 533             if ~discrete
 534                 if or(type(data0)==[1 2 5])
 535                     if grep(normMode, gh+["inf" "all"]+"/", "r")~=[]
 536                         N = N + Ninf
 537                     end
 538                     if grep(normMode, gh+["nan" "all"]+"/", "r")~=[]
 539                         N = N + Nnan
 540                     end
 541                     if type(data0)==5 & ..
 542                        grep(normMode,gh+["zeros" "all"]+"/", "r")~=[]
 543                         N = N + Nzeros
 544                     end
 545                 end
 546                 if grep(normMode, gh+["leftout" "out" "all"]+"/", "r")~=[]
 547                     N = N + Nleftout
 548                 end
 549                 if grep(normMode, gh+["rightout" "out" "all"]+"/", "r")~=[]
 550                     N = N + Nrightout
 551                 end
 552                 if type(data0)==10 & grep(normMode, gh+["empty" "all"]+"/", "r")~=[]
 553                     N = N + Nempty
 554                 end
 555             // Normalized discrete mode:
 556             else
 557                 if grep(normMode, gh+"all/", "r")
 558                     N = size(data0,"*")
 559                 else
 560                     if grep(normMode, gh+"out/", "r")
 561                         N = N + Nout
 562                     end
 563                     if or(type(data0)==[1 2 5])
 564                         if grep(normMode, gh+"nan/", "r")
 565                             N = N + Nnan
 566                         end
 567                         if grep(normMode, gh+"inf/", "r")
 568                             N = N + Ninf
 569                             // Inf already in bins must not be counted twice:
 570                             binsinf = isinf(bins)
 571                             if or(binsinf)
 572                                 N = N - sum(heights(binsinf))
 573                             end
 574                         end
 575                         if type(data0)==5 & and(bins~=0) & ..
 576                            grep(normMode, gh+"zeros/", "r")~=[]
 577                             N = N + Nzeros
 578                         end
 579                     end
 580                     if type(data0)==10 & ~or(bins=="") & ..
 581                         grep(normMode, gh+"empty/", "r")~=[]
 582                         N = N + Nempty
 583                     end
 584                 end
 585             end
 586         //else    // all (but zeros for sparses)
 587         //    N = size(data,"*")
 588         end
 589         heights = heights / N
 590         jokers = jokers / N
 591     end
 592     if grep(histmode, "density") then
 593         heights = heights ./ (bins(2:$) - bins(1:($-1)))
 594     end
 595
 596     // Formatting outputs for sparse input data
 597     // ----------------------------------------
 598     if type(data0)==5 then
 599         ind = sparse(IJ, ind, size(data0))
 600     end
 601
 602     // Shifting argouts in Scilab 5.5 compatibility mode
 603     // -------------------------------------------------
 604     if mode55 then
 605         jokers = ind
 606     end
 607 endfunction