1 /* Algorithm AS 51 Appl. Statist. (1972), vol. 21, p. 218
2 original (C) Royal Statistical Society 1972
3
4 Performs an iterative proportional fit of the marginal totals of a
5 contingency table.
6 */
7
8 #ifdef HAVE_CONFIG_H
9 # include <config.h>
10 #endif
11
12 #include <math.h>
13
14 #include <stdio.h>
15 #include <R_ext/Memory.h>
16 #include <R_ext/Applic.h>
17 #include <R_ext/Error.h> /* for error */
18
19 #undef max
20 #undef min
21 #undef abs
22 #define max(a, b) ((a) < (b) ? (b) : (a))
23 #define min(a, b) ((a) > (b) ? (b) : (a))
24 #define abs(x) ((x) >= 0 ? (x) : -(x))
25
26 static void collap(int nvar, double *x, double *y, int locy,
27 int *dim, int *config);
28 static void adjust(int nvar, double *x, double *y, double *z,
29 int *locz, int *dim, int *config, double *d);
30
31 /* Table of constant values */
32
33 static void
loglin(int nvar,int * dim,int ncon,int * config,int ntab,double * table,double * fit,int * locmar,int nmar,double * marg,int nu,double * u,double maxdev,int maxit,double * dev,int * nlast,int * ifault)34 loglin(int nvar, int *dim, int ncon, int *config, int ntab,
35 double *table, double *fit, int *locmar, int nmar, double *marg,
36 int nu, double *u, double maxdev, int maxit,
37 double *dev, int *nlast, int *ifault)
38 {
39 // nvar could be zero (no-segfault test)
40 if (!nvar) error("no variables"); // not translated
41 int i, j, k, n, point, size, check[nvar], icon[nvar];
42 double x, y, xmax;
43
44 /* Parameter adjustments */
45 --dim;
46 --locmar;
47 config -= nvar + 1;
48 --fit;
49 --table;
50 --marg;
51 --u;
52 --dev;
53
54 /* Function body */
55
56 *ifault = 0;
57 *nlast = 0;
58
59 /* Check validity of NVAR, the number of variables, and of maxit,
60 the maximum number of iterations */
61
62 if (nvar > 0 && maxit > 0) goto L10;
63 L5:
64 *ifault = 4;
65 return;
66
67 /* Look at table and fit constants */
68
69 L10:
70 size = 1;
71 for (j = 1; j <= nvar; j++) {
72 if (dim[j] <= 0) goto L5;
73 size *= dim[j];
74 }
75 if (size <= ntab) goto L40;
76 L35:
77 *ifault = 2;
78 return;
79 L40:
80 x = 0.;
81 y = 0.;
82 for (i = 1; i <= size; i++) {
83 if (table[i] < 0. || fit[i] < 0.) goto L5;
84 x += table[i];
85 y += fit[i];
86 }
87
88 /* Make a preliminary adjustment to obtain the fit to an empty
89 configuration list */
90
91 if (y == 0.) goto L5;
92 x /= y;
93 for (i = 1; i <= size; i++) fit[i] = x * fit[i];
94 if (ncon <= 0 || config[nvar + 1] == 0) return;
95
96 /* Allocate marginal tables */
97
98 point = 1;
99 for (i = 1; i <= ncon; i++) {
100 /* A zero beginning a configuration indicates that the list is
101 completed */
102 if (config[i * nvar + 1] == 0) goto L160;
103 /* Get marginal table size. While doing this task, see if the
104 configuration list contains duplications or elements out of
105 range. */
106 size = 1;
107 for (j = 0; j < nvar; j++) check[j] = 0;
108 for (j = 1; j <= nvar; j++) {
109 k = config[j + i * nvar];
110 /* A zero indicates the end of the string. */
111 if (k == 0) goto L130;
112 /* See if element is valid. */
113 if (k >= 0 && k <= nvar) goto L100;
114 L95:
115 *ifault = 1;
116 return;
117 /* Check for duplication */
118 L100:
119 if (check[k - 1]) goto L95;
120 check[k - 1] = 1;
121 /* Get size */
122 size *= dim[k];
123 }
124
125 /* Since U is used to store fitted marginals, size must not
126 exceed NU */
127 L130:
128 if (size > nu) goto L35;
129
130 /* LOCMAR points to marginal tables to be placed in MARG */
131 locmar[i] = point;
132 point += size;
133 }
134
135 /* Get N, number of valid configurations */
136
137 i = ncon + 1;
138 L160:
139 n = i - 1;
140
141 /* See if MARG can hold all marginal tables */
142
143 if (point > nmar + 1) goto L35;
144
145 /* Obtain marginal tables */
146
147 for (i = 1; i <= n; i++) {
148 for (j = 1; j <= nvar; j++) {
149 icon[j - 1] = config[j + i * nvar];
150 }
151 collap(nvar, &table[1], &marg[1], locmar[i], &dim[1], icon);
152 }
153
154 /* Perform iterations */
155
156 for (k = 1; k <= maxit; k++) {
157 /* XMAX is maximum deviation observed between fitted and true
158 marginal during a cycle */
159 xmax = 0.;
160 for (i = 1; i <= n; i++) {
161 for (j = 1; j <= nvar; j++) icon[j - 1] = config[j + i * nvar];
162 collap(nvar, &fit[1], &u[1], 1, &dim[1], icon);
163 adjust(nvar, &fit[1], &u[1], &marg[1], &locmar[i], &dim[1], icon, &xmax);
164 }
165 /* Test convergence */
166 dev[k] = xmax;
167 if (xmax < maxdev) goto L240;
168 }
169 if (maxit > 1) goto L230;
170 *nlast = 1;
171 return;
172
173 /* No convergence */
174 L230:
175 *ifault = 3;
176 *nlast = maxit;
177 return;
178
179 /* Normal termination */
180 L240:
181 *nlast = k;
182
183 return;
184 }
185
186 /* Algorithm AS 51.1 Appl. Statist. (1972), vol. 21, p. 218
187
188 Computes a marginal table from a complete table.
189 All parameters are assumed valid without test.
190
191 The larger table is X and the smaller one is Y.
192 */
193
collap(int nvar,double * x,double * y,int locy,int * dim,int * config)194 void collap(int nvar, double *x, double *y, int locy, int *dim, int *config)
195 {
196 int i, j, k, l, n, locu, size[nvar + 1], coord[nvar];
197
198 /* Parameter adjustments */
199 --config;
200 --dim;
201 --x;
202 --y;
203
204 /* Initialize arrays */
205
206 size[0] = 1;
207 for (k = 1; k <= nvar; k++) {
208 l = config[k];
209 if (l == 0) goto L20;
210 size[k] = size[k - 1] * dim[l];
211 }
212
213 /* Find number of variables in configuration */
214
215 k = nvar + 1;
216 L20:
217 n = k - 1;
218
219 /* Initialize Y. First cell of marginal table is at Y(LOCY) and
220 table has SIZE(K) elements */
221
222 locu = locy + size[k - 1] - 1;
223 for (j = locy; j <= locu; j++) y[j] = 0.;
224
225 /* Initialize coordinates */
226
227 for (k = 0; k < nvar; k++) coord[k] = 0;
228
229 /* Find locations in tables */
230 i = 1;
231 L60:
232 j = locy;
233 for (k = 1; k <= n; k++) {
234 l = config[k];
235 j += coord[l - 1] * size[k - 1];
236 }
237 y[j] += x[i];
238
239 /* Update coordinates */
240
241 i++;
242 for (k = 1; k <= nvar; k++) {
243 coord[k - 1]++;
244 if (coord[k - 1] < dim[k]) goto L60;
245 coord[k - 1] = 0;
246 }
247
248 return;
249 }
250
251
252 /* Algorithm AS 51.2 Appl. Statist. (1972), vol. 21, p. 218
253
254 Makes proportional adjustment corresponding to CONFIG.
255 All parameters are assumed valid without test.
256 */
257
adjust(int nvar,double * x,double * y,double * z,int * locz,int * dim,int * config,double * d)258 void adjust(int nvar, double *x, double *y, double *z, int *locz,
259 int *dim, int *config, double *d)
260 {
261 int i, j, k, l, n, size[nvar + 1], coord[nvar];
262 double e;
263
264 /* Parameter adjustments */
265 --config;
266 --dim;
267 --x;
268 --y;
269 --z;
270
271 /* Set size array */
272
273 size[0] = 1;
274 for (k = 1; k <= nvar; k++) {
275 l = config[k];
276 if (l == 0) goto L20;
277 size[k] = size[k - 1] * dim[l];
278 }
279
280 /* Find number of variables in configuration */
281
282 k = nvar + 1;
283 L20:
284 n = k - 1;
285
286 /* Test size of deviation */
287
288 l = size[k - 1];
289 j = 1;
290 k = *locz;
291 for (i = 1; i <= l; i++) {
292 e = abs(z[k] - y[j]);
293 if (e > *d) {
294 *d = e;
295 }
296 j++;
297 k++;
298 }
299
300 /* Initialize coordinates */
301
302 for (k = 0; k < nvar; k++) coord[k] = 0;
303 i = 1;
304
305 /* Perform adjustment */
306
307 L50:
308 j = 0;
309 for (k = 1; k <= n; k++) {
310 l = config[k];
311 j += coord[l - 1] * size[k - 1];
312 }
313 k = j + *locz;
314 j++;
315
316 /* Note that Y(J) should be non-negative */
317
318 if (y[j] <= 0.) x[i] = 0.;
319 if (y[j] > 0.) x[i] = x[i] * z[k] / y[j];
320
321 /* Update coordinates */
322
323 i++;
324 for (k = 1; k <= nvar; k++) {
325 coord[k - 1]++;
326 if (coord[k - 1] < dim[k]) goto L50;
327 coord[k - 1] = 0;
328 }
329
330 return;
331 }
332
333 #undef max
334 #undef min
335 #undef abs
336
337 #include <R.h>
338 #include <Rinternals.h>
339 #include "statsR.h"
340 #ifdef ENABLE_NLS
341 #include <libintl.h>
342 #define _(String) dgettext ("stats", String)
343 #else
344 #define _(String) (String)
345 #endif
346
LogLin(SEXP dtab,SEXP conf,SEXP table,SEXP start,SEXP snmar,SEXP eps,SEXP iter)347 SEXP LogLin(SEXP dtab, SEXP conf, SEXP table, SEXP start,
348 SEXP snmar, SEXP eps, SEXP iter)
349 {
350 int nvar = length(dtab),
351 ncon = ncols(conf),
352 ntab = length(table),
353 nmar = asInteger(snmar),
354 maxit = asInteger(iter),
355 nlast, ifault;
356 double maxdev = asReal(eps);
357 SEXP fit = PROTECT(TYPEOF(start) == REALSXP ? duplicate(start) :
358 coerceVector(start, REALSXP)),
359 locmar = PROTECT(allocVector(INTSXP, ncon)),
360 marg = PROTECT(allocVector(REALSXP, nmar)),
361 u = PROTECT(allocVector(REALSXP, ntab)),
362 dev = PROTECT(allocVector(REALSXP, maxit));
363 dtab = PROTECT(coerceVector(dtab, INTSXP));
364 conf = PROTECT(coerceVector(conf, INTSXP));
365 table = PROTECT(coerceVector(table, REALSXP));
366 loglin(nvar, INTEGER(dtab), ncon, INTEGER(conf), ntab,
367 REAL(table), REAL(fit), INTEGER(locmar), nmar, REAL(marg),
368 ntab, REAL(u), maxdev, maxit, REAL(dev), &nlast, &ifault);
369 switch(ifault) {
370 case 1:
371 case 2:
372 error(_("this should not happen")); break;
373 case 3:
374 warning(_("algorithm did not converge")); break;
375 case 4:
376 error(_("incorrect specification of 'table' or 'start'")); break;
377 default:
378 break;
379 }
380
381 SEXP ans = PROTECT(allocVector(VECSXP, 3)), nm;
382 SET_VECTOR_ELT(ans, 0, fit);
383 SET_VECTOR_ELT(ans, 1, dev);
384 SET_VECTOR_ELT(ans, 2, ScalarInteger(nlast));
385 nm = allocVector(STRSXP, 3);
386 setAttrib(ans, R_NamesSymbol, nm);
387 SET_STRING_ELT(nm, 0, mkChar("fit"));
388 SET_STRING_ELT(nm, 1, mkChar("dev"));
389 SET_STRING_ELT(nm, 2, mkChar("nlast"));
390 UNPROTECT(9);
391 return ans;
392 }
393
394