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