| /* Algorithm AS 159 Applied Statistics (1981), vol. 30, no. 1 |
| original (C) Royal Statistical Society 1981 |
| |
| Generate random two-way table with given marginal totals. |
| |
| Heavily pretty edited by Martin Maechler, Dec 2003 |
| use double precision for integer multiplication (against overflow); |
| */ |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| #ifdef ENABLE_NLS |
| #include <libintl.h> |
| #define _(String) dgettext ("stats", String) |
| #else |
| #define _(String) (String) |
| #endif |
| |
| #include <math.h> |
| |
| #include <R_ext/Random.h> |
| #include <R_ext/Applic.h> |
| #include <R_ext/Boolean.h> |
| #include <R_ext/Error.h> |
| #include <R_ext/Utils.h> |
| |
| void |
| rcont2(int *nrow, int *ncol, |
| /* vectors of row and column totals, and their sum ntotal: */ |
| int *nrowt, int *ncolt, int *ntotal, |
| double *fact, int *jwork, int *matrix) |
| { |
| int j, l, m, ia, ib, ic, jc, id, ie, ii, nll, nlm, nr_1, nc_1; |
| double x, y, dummy, sumprb; |
| Rboolean lsm, lsp; |
| |
| nr_1 = *nrow - 1; |
| nc_1 = *ncol - 1; |
| |
| ib = 0; /* -Wall */ |
| |
| /* Construct random matrix */ |
| for (j = 0; j < nc_1; ++j) |
| jwork[j] = ncolt[j]; |
| |
| jc = *ntotal; |
| |
| for (l = 0; l < nr_1; ++l) { /* ----- matrix[ l, * ] ----- */ |
| ia = nrowt[l]; |
| ic = jc; |
| jc -= ia;/* = n_tot - sum(nr[0:l]) */ |
| |
| for (m = 0; m < nc_1; ++m) { |
| id = jwork[m]; |
| ie = ic; |
| ic -= id; |
| ib = ie - ia; |
| ii = ib - id; |
| |
| if (ie == 0) { /* Row [l,] is full, fill rest with zero entries */ |
| for (j = m; j < nc_1; ++j) |
| matrix[l + j * *nrow] = 0; |
| ia = 0; |
| break; |
| } |
| |
| /* Generate pseudo-random number */ |
| dummy = unif_rand(); |
| |
| do {/* Outer Loop */ |
| |
| /* Compute conditional expected value of MATRIX(L, M) */ |
| |
| nlm = (int)(ia * (id / (double) ie) + 0.5); |
| x = exp(fact[ia] + fact[ib] + fact[ic] + fact[id] |
| - fact[ie] - fact[nlm] |
| - fact[id - nlm] - fact[ia - nlm] - fact[ii + nlm]); |
| if (x >= dummy) |
| break; |
| if (x == 0.)/* MM: I haven't seen this anymore */ |
| error(_("rcont2 [%d,%d]: exp underflow to 0; algorithm failure"), l, m); |
| |
| sumprb = x; |
| y = x; |
| nll = nlm; |
| |
| do { |
| /* Increment entry in row L, column M */ |
| j = (int)((id - nlm) * (double)(ia - nlm)); |
| lsp = (j == 0); |
| if (!lsp) { |
| ++nlm; |
| x = x * j / ((double) nlm * (ii + nlm)); |
| sumprb += x; |
| if (sumprb >= dummy) |
| goto L160; |
| } |
| |
| do { |
| R_CheckUserInterrupt(); |
| |
| /* Decrement entry in row L, column M */ |
| j = (int)(nll * (double)(ii + nll)); |
| lsm = (j == 0); |
| if (!lsm) { |
| --nll; |
| y = y * j / ((double) (id - nll) * (ia - nll)); |
| sumprb += y; |
| if (sumprb >= dummy) { |
| nlm = nll; |
| goto L160; |
| } |
| /* else */ |
| if (!lsp) |
| break;/* to while (!lsp) */ |
| } |
| } while (!lsm); |
| |
| } while (!lsp); |
| |
| dummy = sumprb * unif_rand(); |
| |
| } while (1); |
| |
| L160: |
| matrix[l + m * *nrow] = nlm; |
| ia -= nlm; |
| jwork[m] -= nlm; |
| } |
| matrix[l + nc_1 * *nrow] = ia;/* last column in row l */ |
| } |
| |
| /* Compute entries in last row of MATRIX */ |
| for (m = 0; m < nc_1; ++m) |
| matrix[nr_1 + m * *nrow] = jwork[m]; |
| |
| matrix[nr_1 + nc_1 * *nrow] = ib - matrix[nr_1 + (nc_1-1) * *nrow]; |
| |
| return; |
| } |