The R Project SVN R

Rev

Rev 24068 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 24068 Rev 27560
Line 1... Line 1...
1 
/* Algorithm AS 159 Applied Statistics (1981), vol. 30, no. 1
 1 
/* Algorithm AS 159 Applied Statistics (1981), vol. 30, no. 1
2 
   original (C) Royal Statistical Society 1981
 2 
   original (C) Royal Statistical Society 1981
3 
 
 3 
 
4 
   Generate random two-way table with given marginal totals.
 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);
5 
   */
 8 
*/
6 
 
 9 
 
7 
#include <math.h>
 10 
#include <math.h>
8 
 
 11 
 
9 
#include <R_ext/Random.h>
 12 
#include <R_ext/Random.h>
- 
 
 13 
#include <R_ext/Applic.h>
- 
 
 14 
#include <R_ext/Boolean.h>
- 
 
 15 
#include <R_ext/Error.h>
- 
 
 16 
#include <R_ext/Utils.h>
10 
 
 17 
 
11 
void
 18 
void
- 
 
 19 
rcont2(int *nrow, int *ncol,
- 
 
 20 
       /* vectors of row and column totals, and their sum ntotal: */
12 
rcont2(int *nrow, int *ncol, int *nrowt, int *ncolt, int *ntotal,
 21 
       int *nrowt, int *ncolt, int *ntotal,
13 
       double *fact, int *jwork, int *matrix)
 22 
       double *fact, int *jwork, int *matrix)
14 
{
 23 
{
15 
    int nlmp, j, l, m, ia, ib, ic, jc, id, ie, ii, nrowtl, iap, idp,
 24 
    int j, l, m, ia, ib, ic, jc, id, ie, ii, nll, nlm, nr_1, nc_1;
16 
	igp, ihp, iip, nll, nlm, nrowm, ncolm, lsm, lsp;
 - 
 
17 
    double x, y, dummy, sumprb;
 25 
    double x, y, dummy, sumprb;
- 
 
 26 
    Rboolean lsm, lsp;
18 
 
 27 
 
19 
    --jwork;
 - 
 
20 
    --ncolt;
 - 
 
21 
    --nrowt;
 - 
 
22 
    matrix -= *nrow + 1;
 - 
 
23 
 
 - 
 
24 
    nrowm = *nrow - 1;
 28 
    nr_1 = *nrow - 1;
25 
    ncolm = *ncol - 1;
 29 
    nc_1 = *ncol - 1;
26 
 
 30 
 
27 
    ib = 0;			/* -Wall */
 31 
    ib = 0; /* -Wall */
28 
 
 32 
 
29 
    /* Construct random matrix */
 33 
    /* Construct random matrix */
30 
    for (j = 1; j <= ncolm; ++j) {
 34 
    for (j = 0; j < nc_1; ++j)
31 
	jwork[j] = ncolt[j];
 35 
	jwork[j] = ncolt[j];
32 
    }
 36 
 
33 
    jc = *ntotal;
 37 
    jc = *ntotal;
34 
 
 38 
 
35 
    /* HOP = ONE */
 - 
 
36 
    for (l = 1; l <= nrowm; ++l) {
 39 
    for (l = 0; l < nr_1; ++l) { /* -----  matrix[ l, * ] ----- */
37 
	nrowtl = nrowt[l];
 - 
 
38 
	ia = nrowtl;
 40 
	ia = nrowt[l];
39 
	ic = jc;
 41 
	ic = jc;
40 
	jc -= nrowtl;
 42 
	jc -= ia;/* = n_tot - sum(nr[0:l]) */
- 
 
 43 
 
41 
	for (m = 1; m <= ncolm; ++m) {
 44 
	for (m = 0; m < nc_1; ++m) {
42 
	    id = jwork[m];
 45 
	    id = jwork[m];
43 
	    ie = ic;
 46 
	    ie = ic;
44 
	    ic -= id;
 47 
	    ic -= id;
45 
	    ib = ie - ia;
 48 
	    ib = ie - ia;
46 
	    ii = ib - id;
 49 
	    ii = ib - id;
47 
 
 50 
 
48 
	    /* Test for zero entries in MATRIX */
 51 
	    if (ie == 0) { /* Row [l,] is full, fill rest with zero entries */
49 
	    if (ie != 0) {
 52 
		for (j = m; j < nc_1; ++j)
- 
 
 53 
		    matrix[l + j * *nrow] = 0;
50 
		goto L130;
 54 
		ia = 0;
- 
 
 55 
		break;
51 
	    }
 56 
	    }
52 
	    for (j = m; j <= *ncol; ++j) {
 - 
 
53 
		matrix[l + j * *nrow] = 0;
 - 
 
54 
	    }
 - 
 
55 
	    goto L190;
 - 
 
56 
 
 57 
 
57 
	    /* Generate pseudo-random number */
 58 
	    /* Generate pseudo-random number */
58 
L130:
 - 
 
59 
	    dummy = unif_rand();
 59 
	    dummy = unif_rand();
60 
 
 60 
 
- 
 
 61 
	    do {/* Outer Loop */
- 
 
 62 
 
61 
	    /* Compute conditional expected value of MATRIX(L, M) */
 63 
		/* Compute conditional expected value of MATRIX(L, M) */
62 
L131:
 64 
 
63 
	    nlm = ia * id / (double) ie + 0.5;
 65 
		nlm = ia * (id / (double) ie) + 0.5;
- 
 
 66 
		x = exp(fact[ia] + fact[ib] + fact[ic] + fact[id]
64 
	    iap = ia + 1;
 67 
			- fact[ie] - fact[nlm]
- 
 
 68 
			- fact[id - nlm] - fact[ia - nlm] - fact[ii + nlm]);
65 
	    idp = id + 1;
 69 
		if (x >= dummy)
66 
	    igp = idp - nlm;
 70 
		    break;
- 
 
 71 
		if (x == 0.)/* MM: I haven't seen this anymore */
- 
 
 72 
		    error("rcont2 [%d,%d]: exp underflow to 0;"
67 
	    ihp = iap - nlm;
 73 
			  " algorithm failure", l, m);
- 
 
 74 
 
68 
	    nlmp = nlm + 1;
 75 
		sumprb = x;
- 
 
 76 
		y = x;
69 
	    iip = ii + nlmp;
 77 
		nll = nlm;
- 
 
 78 
 
- 
 
 79 
		do {
70 
	    x = exp(fact[iap - 1] + fact[ib] + fact[ic] +
 80 
		    /* Increment entry in row L, column M */
71 
		    fact[idp - 1] - fact[ie] - fact[nlmp - 1] -
 81 
		    j = (id - nlm) * (double)(ia - nlm);
72 
		    fact[igp - 1] - fact[ihp - 1] - fact[iip - 1]);
 82 
		    lsp = (j == 0);
73 
	    if (x >= dummy) {
 83 
		    if (!lsp) {
74 
		goto L160;
 84 
			++nlm;
- 
 
 85 
			x = x * j / ((double) nlm * (ii + nlm));
75 
	    }
 86 
			sumprb += x;
76 
	    sumprb = x;
 87 
			if (sumprb >= dummy)
77 
	    y = x;
 88 
			    goto L160;
78 
	    nll = nlm;
 89 
		    }
- 
 
 90 
 
79 
	    lsp = (0);
 91 
		    do {
80 
	    lsm = (0);
 92 
			R_CheckUserInterrupt();
81 
 
 93 
 
82 
	    /* Increment entry in row L, column M */
 94 
			/* Decrement entry in row L, column M */
83 
L140:
 - 
 
84 
	    j = (id - nlm) * (ia - nlm);
 95 
			j = nll * (double)(ii + nll);
85 
	    if (j == 0) {
 96 
			lsm = (j == 0);
86 
		goto L156;
 97 
			if (!lsm) {
87 
	    }
 - 
 
88 
	    ++nlm;
 98 
			    --nll;
89 
	    x = x * j / (double) (nlm * (ii + nlm));
 99 
			    y = y * j / ((double) (id - nll) * (ia - nll));
90 
	    sumprb += x;
 100 
			    sumprb += y;
91 
	    if (sumprb >= dummy) {
 101 
			    if (sumprb >= dummy) {
- 
 
 102 
				nlm = nll;
92 
		goto L160;
 103 
				goto L160;
93 
	    }
 104 
			    }
- 
 
 105 
			    /* else */
- 
 
 106 
			    if (!lsp)
- 
 
 107 
				break;/* to while (!lsp) */
94 
L150:
 108 
			}
95 
	    if (lsm) {
 109 
		    } while (!lsm);
- 
 
 110 
 
96 
		goto L155;
 111 
		} while (!lsp);
- 
 
 112 
 
- 
 
 113 
		dummy = sumprb * unif_rand();
- 
 
 114 
 
97 
	    }
 115 
	    } while (1);
98 
 
 116 
 
99 
	    /* Decrement entry in row L, column M */
 - 
 
100 
	    j = nll * (ii + nll);
 - 
 
101 
	    if (j == 0) {
 - 
 
102 
		goto L154;
 - 
 
103 
	    }
 - 
 
104 
	    --nll;
 - 
 
105 
	    y = y * j / (double) ((id - nll) * (ia - nll));
 - 
 
106 
	    sumprb += y;
 - 
 
107 
	    if (sumprb >= dummy) {
 - 
 
108 
		goto L159;
 - 
 
109 
	    }
 - 
 
110 
	    if (! lsp) {
 - 
 
111 
		goto L140;
 - 
 
112 
	    }
 - 
 
113 
	    goto L150;
 - 
 
114 
L154:
 - 
 
115 
	    lsm = (1);
 - 
 
116 
L155:
 - 
 
117 
	    if (!lsp) {
 - 
 
118 
		goto L140;
 - 
 
119 
	    }
 - 
 
120 
	    dummy = sumprb * unif_rand();
 - 
 
121 
	    goto L131;
 - 
 
122 
L156:
 - 
 
123 
	    lsp = (1);
 - 
 
124 
	    goto L150;
 - 
 
125 
L159:
 - 
 
126 
	    nlm = nll;
 - 
 
127 
L160:
 117 
L160:
128 
	    matrix[l + m * *nrow] = nlm;
 118 
	    matrix[l + m * *nrow] = nlm;
129 
	    ia -= nlm;
 119 
	    ia -= nlm;
130 
	    jwork[m] -= nlm;
 120 
	    jwork[m] -= nlm;
131 
	}
 121 
	}
132 
	matrix[l + *ncol * *nrow] = ia;
 122 
	matrix[l + nc_1 * *nrow] = ia;/* last column in row l */
133 
L190:
 - 
 
134 
	;
 - 
 
135 
    }
 123 
    }
136 
 
 124 
 
137 
    /* Compute entries in last row of MATRIX */
 125 
    /* Compute entries in last row of MATRIX */
- 
 
 126 
    for (m = 0; m < nc_1; ++m)
- 
 
 127 
	matrix[nr_1 + m * *nrow] = jwork[m];
138 
 
 128 
 
139 
    for (m = 1; m <= ncolm; ++m) {
 - 
 
140 
	matrix[*nrow + m * *nrow] = jwork[m];
 - 
 
141 
    }
 - 
 
142 
    matrix[*nrow + *ncol * *nrow] = ib - matrix[*nrow + ncolm * *nrow];
 129 
    matrix[nr_1 + nc_1 * *nrow] = ib - matrix[nr_1 + (nc_1-1) * *nrow];
143 
 
 130 
 
144 
    return;
 131 
    return;
145 
}
 132 
}