1 #include <math.h>
2 #include <R.h>
3
4 #define max(A,B) ((A) > (B) ? (A):(B))
5 #define min(A,B) ((A) < (B) ? (A):(B))
6
GMLE(int * Mstrata,int * Istrata,int * Mindex,int * Iindex,int * N,int * M,double * z,double * oldZ,double * tol,int * maxstep,int * niter)7 void GMLE(int *Mstrata,
8 int *Istrata,
9 int *Mindex,
10 int *Iindex,
11 int *N,
12 int *M,
13 double *z,
14 double *oldZ,
15 double *tol,
16 int *maxstep,
17 int *niter){
18
19 int i,j,k,l,m,step,done;
20 double newZ,nom, denom, diff;
21 step=0;
22 done=0;
23 while (done==0 && step < *maxstep){
24 /* Rprintf("\n\nStep=%d\t\n",step); */
25 diff=0;
26 for(k=0;k<*M;k++) oldZ[k]= z[k];
27 for(k=0;k<*M;k++){
28 nom=0;
29 newZ=0;
30 for(j=Mstrata[k]; j< Mstrata[k+1];j++){
31 i=Mindex[j]-1;
32 denom=0;
33 for(l=Istrata[i]; l < Istrata[i+1];l++){
34 m=Iindex[l]-1;
35 denom += oldZ[m];
36 }
37 nom = oldZ[k];
38 newZ += nom/denom;
39 }
40 z[k]=newZ/(*N);
41 }
42 for (k=0;k<*M;k++){
43 /* Rprintf("k=%d\toldZ[k]=%1.2f\tz[k]=%1.2f\tdiff=%1.2f\t\n",k,oldZ[k],z[k],diff); */
44 diff=max(max(z[k]-oldZ[k],oldZ[k]-z[k]),diff);
45 }
46 if (diff < *tol) done=1;
47 step++;
48 }
49 niter[0]=step;
50 }
51