trandcounts: add repeat option - numtools - perform numerical operations on vectors and matrices in unix pipes
 (HTM) git clone git://src.adamsgaard.dk/numtools
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit 9dc61cb5c221a9fef82704d7edd23c8f61608ca1
 (DIR) parent 274c96b0785a902cc65d91cdf191593cdac7a4f5
 (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
       Date:   Mon,  9 May 2022 10:09:02 +0200
       
       randcounts: add repeat option
       
       Diffstat:
         M randcounts.1                        |      12 +++++++++---
         M randcounts.c                        |      43 +++++++++++++++++--------------
       
       2 files changed, 33 insertions(+), 22 deletions(-)
       ---
 (DIR) diff --git a/randcounts.1 b/randcounts.1
       t@@ -8,6 +8,7 @@
        .Nm
        .Op Fl h
        .Op Fl n Ar num
       +.Op Fl r Ar repeats
        .Op Fl s Ar seed
        .Ar weight1
        .Op Ar weight2 ...
       t@@ -34,6 +35,8 @@ Show usage information.
        .It Fl n Ar num
        Number of points to place in the bins.
        The default is 1.
       +.It Fl r Ar repeats
       +Repeat the binning several times, with one realization per line of output.
        .It Fl s Ar seed
        Seed the pseudo-random number generator with this value, which is used
        to generate reproducable binning.
       t@@ -48,11 +51,14 @@ Due to the randomness, your output may differ:
        Put 100 points in two bins with 75% and 25% probability, respectively:
        .Pp
        .Dl $ randcounts -n 100 0.75 0.25
       -.Dl 78        22
       +.Pp
       +Put 100 points in three equal bins 1000 times, and calculate the average bin sizes with
       +.Xr mean 1 :
       +.Pp
       +.Dl $ randcounts -n 100 -r 1000 0.333 0.333 0.334 | mean
       +.Dl 33.067        32.82        34.113
        .Sh SEE ALSO
       -.Xr max 1 ,
        .Xr mean 1 ,
       -.Xr min 1 ,
        .Xr range 1 ,
        .Xr rangetest 1 ,
        .Xr sum 1
 (DIR) diff --git a/randcounts.c b/randcounts.c
       t@@ -14,7 +14,7 @@ char *argv0;
        static void
        usage(void)
        {
       -        errx(1, "usage: %s [-h] [-n num] [-s seed] "
       +        errx(1, "usage: %s [-h] [-n num] [-r repeats] [-s seed] "
                        "weight1 [weight2 ...]\n", argv0);
        }
        
       t@@ -22,7 +22,7 @@ int
        main(int argc, char *argv[])
        {
                int i, ret, s = 0, N;
       -        long j, seed, *counts = NULL, n = 1;
       +        long j, seed, *counts = NULL, n = 1, r, repeats = 1;
                double val = 0.0, weightsum = 0.0, *weights = NULL;
                struct timeval t1;
        
       t@@ -36,6 +36,9 @@ main(int argc, char *argv[])
                case 'n':
                        n = atol(EARGF(usage()));
                        break;
       +        case 'r':
       +                repeats = atol(EARGF(usage()));
       +                break;
                case 's':
                        s = 1;
                        seed = atol(EARGF(usage()));
       t@@ -53,7 +56,6 @@ main(int argc, char *argv[])
                        err(1, "calloc");
        
                for (i = 0; i < N; i++) {
       -                counts[i] = 0;
                        weights[i] = atof(argv[i]);
                        if (weights[i] <= 0.0)
                                errx(1, "weight %d is not positive (%g)", i, weights[i]);
       t@@ -77,24 +79,27 @@ main(int argc, char *argv[])
                }
        #endif
        
       -        for (j = 0; j < n; j++) {
       -                val = drand48();
       -                weightsum = 0.0;
       -                for (i = 0; i < N; i++) {
       -                        weightsum += weights[i];
       -                        if (val <= weightsum) {
       -                                counts[i]++;
       -                                break;
       +        for (r = 0; r < repeats; r++) {
       +                for (i = 0; i < N; i++)
       +                        counts[i] = 0;
       +                for (j = 0; j < n; j++) {
       +                        val = drand48();
       +                        weightsum = 0.0;
       +                        for (i = 0; i < N; i++) {
       +                                weightsum += weights[i];
       +                                if (val <= weightsum) {
       +                                        counts[i]++;
       +                                        break;
       +                                }
                                }
                        }
       -        }
       -
       -        for (i = 0; i < N; i++) {
       -                printf("%ld", counts[i]);
       -                if (i < N - 1)
       -                        putchar('\t');
       -                else
       -                        putchar('\n');
       +                for (i = 0; i < N; i++) {
       +                        printf("%ld", counts[i]);
       +                        if (i < N - 1)
       +                                putchar('\t');
       +                        else
       +                                putchar('\n');
       +                }
                }
        
                free(counts);