1 /* Algorithm AS 159 Applied Statistics (1981), vol. 30, no. 1
2 original (C) Royal Statistical Society 1981
3
4 Generate random two-way table with given marginal totals.
5
6 Heavily pretty edited by Martin Maechler, Dec 2003
7 use double precision for integer multiplication (against overflow);
8 */
9
10 #ifdef HAVE_CONFIG_H
11 #include <config.h>
12 #endif
13
14 #ifdef ENABLE_NLS
15 #include <libintl.h>
16 #define _(String) dgettext ("stats", String)
17 #else
18 #define _(String) (String)
19 #endif
20
21 #include <math.h>
22
23 #include <R_ext/Random.h>
24 #include <R_ext/Applic.h>
25 #include <R_ext/Boolean.h>
26 #include <R_ext/Error.h>
27 #include <R_ext/Print.h>
28 #include <R_ext/Utils.h>
29 #ifdef DEBUG_rcont2
30 # include <limits.h>
31 #endif
32
33 #include "stats.h"
34
35 // NB: Exported via S_rcont() --> ../../../include/R_ext/stats_stubs.h & stats_package.h
36 void
rcont2(int nrow,int ncol,const int nrowt[],const int ncolt[],int ntotal,const double fact[],int * jwork,int * matrix)37 rcont2(int nrow, int ncol,
38 /* vectors of row and column totals, and their sum ntotal: */
39 const int nrowt[], const int ncolt[], int ntotal,
40 const double fact[],
41 int *jwork, int *matrix)
42 {
43 int nr_1 = nrow - 1,
44 nc_1 = ncol - 1,
45 ib = 0; /* -Wall */
46
47 /* Construct random matrix */
48 for (int j = 0; j < nc_1; ++j)
49 jwork[j] = ncolt[j];
50
51 int jc = ntotal;
52 for (int l = 0; l < nr_1; ++l) { /* ----- matrix[ l, * ] ----- */
53 int ia = nrowt[l],
54 ic = jc;
55 jc -= ia;/* = n_tot - sum(nr[0:l]) */
56
57 for (int m = 0; m < nc_1; ++m) {
58 int id = jwork[m],
59 ie = ic, ii;
60 ib = ie - ia;
61 ii = ib - id;
62 ic -= id;
63
64 if (ie == 0) { /* Row [l,] is full, fill rest with zero entries */
65 for (int j = m; j < nc_1; ++j)
66 matrix[l + j * nrow] = 0;
67 ia = 0;
68 break;
69 }
70
71 /* Generate pseudo-random number */
72 double U = unif_rand();
73 int nlm;
74 do {/* Outer Loop */
75
76 /* Compute conditional expected value of MATRIX(L, M) */
77
78 nlm = (int)(ia * (id / (double) ie) + 0.5);
79 double x = exp(fact[ia] + fact[ib] + fact[ic] + fact[id]
80 - fact[ie] - fact[nlm]
81 - fact[id - nlm] - fact[ia - nlm] - fact[ii + nlm]);
82 if (x >= U)
83 break;
84 if (x == 0.)/* MM: I haven't seen this anymore */
85 error(_("rcont2 [%d,%d]: exp underflow to 0; algorithm failure"), l, m);
86
87 double sumprb = x,
88 y = x;
89
90 int nll = nlm;
91 Rboolean lsp;
92 do {
93 /* Increment entry in row L, column M */
94 double j = (id - nlm) * (double)(ia - nlm);
95 #ifdef DEBUG_rcont2
96 if(j > INT_MAX)
97 REprintf("Incr.: j = %20.20g > INT_MAX !! (id,ia,nlm) = (%d,%d,%d)\n",
98 j, id,ia,nlm);
99 #endif
100 lsp = ((nlm == ia) || (nlm == id));
101 if (!lsp) {
102 ++nlm;
103 x *= j / ((double) nlm * (ii + nlm));
104 sumprb += x;
105 if (sumprb >= U)
106 goto L160;
107 }
108
109 Rboolean lsm;
110 do {
111 R_CheckUserInterrupt();
112
113 /* Decrement entry in row L, column M */
114 j = (nll * (double)(ii + nll));
115 #ifdef DEBUG_rcont2
116 if(j > INT_MAX)
117 REprintf("Decr.: j = %20.20g > INT_MAX !! (ii,nll) = (%d,%d)\n",
118 j, ii,nll);
119 #endif
120 lsm = (nll == 0);
121 if (!lsm) {
122 --nll;
123 y *= j / ((double) (id - nll) * (ia - nll));
124 sumprb += y;
125 if (sumprb >= U) {
126 nlm = nll;
127 goto L160;
128 }
129 /* else */
130 if (!lsp)
131 break;/* to while (!lsp) */
132 }
133 } while (!lsm);
134
135 } while (!lsp);
136
137 U = sumprb * unif_rand();
138
139 } while (1); // 'Outer Loop'
140
141 L160:
142 matrix[l + m * nrow] = nlm;
143 ia -= nlm;
144 jwork[m] -= nlm;
145 }// for (m = 0..nc_1-1)
146 matrix[l + nc_1 * nrow] = ia;/* last column in row l */
147 } // for (l = ...)
148
149 /* Compute entries in last row of MATRIX */
150 for (int m = 0; m < nc_1; ++m)
151 matrix[nr_1 + m * nrow] = jwork[m];
152
153 matrix[nr_1 + nc_1 * nrow] = ib - matrix[nr_1 + (nc_1-1) * nrow];
154
155 return;
156 }
157