#undef M #undef N #define max(x,y) (x > y) ? x : y #define min(x,y) (x < y) ? x : y #define sq(x) x*x #include #include #include "rng_afm.c" /**************************************************************************** ct.chisq.c current version 8/15/11 For OPTIONC = NO, ct.chisq.c estimates the number of KR X KC contingency tables with given row and column sums and the number of tables with chi-squared statistics <= TEST_STAT, where TEST_STAT = the initial table's chi-squared statistic. The program computes TEST_STAT. For OPTIONC = YES, ct.chisq.c uses an importance sampling techniqu to estimate the number of tables with chi-squared statistics <= TEST_STAT. The program is a modification of ct.hm.c and is based on the multistage MCMC sampling plan described in Fishman, G. S. (2011). Countng subsets of contingency tables, Deparment of Statistics and Operations Research, University of North Carolina, Chapel Hill. For given row and colum sums, it offers two alternative input modes. 1. Built-in instance of a table with the given constraints (see examples below). 2. Scanning a file called ct.hm.initial.txt containing a table satisfying the constraints. It also requires numerical values for the variables: READ = NO if mode 1 is selected = YES if mode 2 is selected KR = no. of rows KC = no. of columns NREPS = no. of independent replications T = sample-path length WU = no. of iterations in warmup interval SW = stage on switch from CTLARGE to CTSMALL occurs OPTIONC = as above if OPTIONC = YES, the program requires values for KSTAR and calls ct_g_max.c to compute GGMAX := hueristic lower bound for g_max(A) as described in the referecne. See the code section below entitled "Assigning parametric values." It also contains other variables with preassigned values that have worked successfully in the nine examples in the technical report. ct.chisq.c also requires rng_amf.c and seed_file. rng_afm.c generates pseudorandom numbers from the Mersenne Twister with initial seeds in seed_file. Before terminating execution, ct.hm.c writes the last table that satisfies the constraints in ct.hm.initial.txt and the last pseudorandom numbers in seed_file. The file ct.hm.out.txt contains the output. *************************************************************************************************/ int main(void) { /* begin ct.chisq.c procedure */ double ct_g_max(int KR, int KC, double MR[], double MC[]); double NSTAGES,DUM3,DUM4,BZ,C2SEZ,C2REZ,EE,ES2A,E2REL,E2DOT,E2REZ,E2SEZ; double CSEZ,CREL,CREZ,DS,CVAR,CS[101],CHIZ,CUM,AZ,PPP,ZZ,FNUM,BIGQ,RR,Z1,DUM1,DUM2; double DE,DS2A,C2REL,C2DOT,AC,CRATIO,ND,CMIN[301],ZF[301]; double SUMXI[301],XI,RHO[301],SA,SG,S2A,DA; double VAR,REL,SEZ,REZ,SUMB,LAM,A[301][1001],HMR[301]; double DOT,DUM,H,HNEW,HOLD,MU,Z,ZBAR,SE,RE; double U,SUM,FI,SLOG,DDD,SSL; double COLD,STD[101][101],STN[101][101],SN,SD,CHLOG,SUC,SUD,CHZ,Z2; double B[301],BR,THETA,PSI; double TEST_STAT,GAM2,GAM,Y[51][51],SC[51],SR[51]; double RMU[301],CMU[301],NUM; double INCX,INCY,CNEW,AX,AY,BX,BY,LX,LY,UX,UY,RADX,RADY; double HMAX,GMAX,T,STEP,OLDC; double SUMXIA[301],XIA,LAMA; int READ,YES=1,NO=0,OPTIONC,W[301],DELTA,QQQ,PICK; int KK,SW,KC,KR,MM,NN,KZ[101],R; int I,J,K,IP,JP,RI,RJ,NREPS; int LSTAGE,STAGE,KSTAR,L,LL,WU; /* A.1 Diaconis and Holmes double GGMAX = 253.5785; double X[6][4] = {{0,0,0,0}, {0,5,2,3}, {0,50,7,5}, {0,3,6,4}, {0,5,3,3}, {0,2,7,30}}; */ /* A.2 Koch et a. (1983). Effectiveness of analgesic drug double GGMAX = 357.3234; double X[9][4] = {{0,0,0,0}, {0,3,20,5}, {0,11,14,8}, {0,3,14,12}, {0,6,13,5}, {0,12,12,0}, {0,11,10,0}, {0,3,9,4}, {0,6,9,3}}; */ /* A.3 Eye-hair data double GGMAX = 1306.1; double X[5][5] = {{0,0,0,0,0}, {0,68,119,26,7}, {0,20,84,17,94}, {0,15,54,14,10}, {0,5,29,14,16}}; */ /* A.4a SAS eye hair example double GGMAX = 1495.9; double X[4][6] = {{0,0,0,0,0,0}, {0,69,28,68,51, 6}, {0,69,38,55,37, 0}, {0,90,47,94,94,16}}; */ /* A.4b SAS eye hair example double GGMAX = 1495.9; double X[6][4] ={{0,0,0,0}, {0,69,69,90}, {0,28,38,47}, {0,68,55,94}, {0,51,37,94}, {0,6,0,16}}; */ /* A.5 Fienberg table 2-6a KR = 4, KC = 4 p. 20 double GGMAX = 5386.0; int X[5][5] = {{0,0,0,0,0}, {0,239,309,233,53}, {0,6,11,70,199}, {0,1,7,12,215}, {0,794,781,922,501}}; */ /* A.6 Fienberg table 2-6c KR = 5, KC = 4 p. 20p. double GGMAX = 10446; int X[6][5] = {{0,0,0,0,0}, {0,215,208,138,83}, {0,281,285,284,197}, {0,372,386,446,385}, {0,128,176,238,186}, {0,44,53,131,117}}; */ /* A7 Diaconis and Efron double GGMAX =48143; double X[6][5]={{0,0,0,0,0}, {0,2161,3577,2184,1636}, {0,2755,5081,2222,1052}, {0,936,1753,640,306}, {0,225,419,96,38}, {0,39,98,31,14}}; */ /* A8 Fienberg table 2-8 KR = 8, KC = 8 p. 24 */ double GGMAX = 391270; double X[9][9] = {{0,0,0,0,0,0,0,0,0}, {0,26,50,11,6,82,39,48,11}, {0,65,2997,238,85,2553,1083,1349,216}, {0,12,279,197,36,459,197,221,47}, {0,3,102,40,61,243,115,101,38}, {0,75,2628,413,229,12137,2658,3689,687}, {0,52,1117,191,102,2649,3210,1973,301}, {0,42,1251,206,117,3757,1962,4646,391}, {0,3,221,51,24,678,301,367,269}}; /* A.9 Table in Andrews and Herzberg (1985, p. 432) double GGMAX = 798.5538; double X[13][13] = {{0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,1,0,0,0,1,2,0,0,1,0,1,0}, {0,1,0,0,1,0,0,0,0,0,1,0,2}, {0,1,0,0,0,2,1,0,0,0,0,0,1}, {0,3,0,2,0,0,0,1,0,1,3,1,1}, {0,2,1,1,1,1,1,1,1,1,1,1,0}, {0,2,0,0,0,1,0,0,0,0,0,0,0}, {0,2,0,2,1,0,0,0,0,1,1,1,2}, {0,0,0,0,3,0,0,1,0,0,1,0,2}, {0,0,0,0,1,1,0,0,0,0,0,1,0}, {0,1,1,0,2,0,0,1,0,0,1,1,0}, {0,0,1,1,1,2,0,0,2,0,1,1,0}, {0,0,1,1,0,0,0,1,0,0,0,0,0}}; */ /* A.10 Table in Andrews and Herzberg (1985, p. 431) double GGMAX = 4.097437e+03; double X[13][13] = {{0,0,0,0,0,0,0,0,0,0,0,0,0}, {0,2,6,7,4,3,4,1,7,1,2,3,3}, {0,3,4,1,3,4,3,1,1,2,5,1,6}, {0,2,4,7,3,3,1,4,4,1,1,4,3}, {0,0,3,0,2,4,2,2,5,0,2,1,2}, {0,1,3,3,2,1,2,1,5,5,2,1,1}, {0,4,1,0,2,7,1,1,2,1,3,4,0}, {0,5,2,5,1,2,5,2,3,1,2,0,3}, {0,4,3,4,3,4,3,1,4,2,2,3,5}, {0,1,4,4,3,1,0,3,3,7,0,3,2}, {0,2,5,7,2,1,2,5,1,1,3,2,2}, {0,3,6,2,4,4,2,4,3,1,1,5,2}, {0,1,5,3,4,6,2,3,2,5,2,5,2}}; */ FILE *fp,*fp1,*fpa,*fpaa,*fpab; /**************Load seeds for the mersenne twister*************************/ fload_seeds(); /**************************************************************************/ fp = fopen("ct.chisq.out.txt", "w"); fp1 = fopen("ct.hm.initial.txt", "r"); /*************************************************************************** Assign parametric values ***************************************************************************/ READ = NO; /* if READ = YES, read new table from disk. If READ = NO, use above table */ KR = 8; /* no. of rows */ KC = 8; /* no. of columns */ NREPS = 10; /* no. of independent replications */ T = 100000; /* sample-path length */ WU = 1000; /* warm-up interval */ SW = 36; /* if STAGE < SW, use CTLARGE; otherwise, use CTSMALL */ LSTAGE = 1; /* sum of all row sums = sum of all column sums */ NSTAGES = 40; /* no. of stages for simulated annealing */; PSI = 1; /* upper bound constraint on XI_1 and XI_R */ THETA = 1.0; /* growth rate for geometrically increasing increments */ OPTIONC = NO; /* if OPTIONC = YES, estimates # <= TEST_STATISTIC by method C; if OPTIONC = NO, estimates # <= TEST_STATISTIC by methods A and B and total count */ KSTAR = 40; /* stage on which for KSTAR =0, GAM = TEST_STAT. For KSTAR no., scaling on stages 1 through NSTAGES -1 and GAM = TEST_STAT thereafter */ /* If OPTIONC = YES, then GGMAX also needs a value. For ten examples, the corresponding values are listed above. For a new problem it may be calculated using the ~gmax.c C program. */ /********************************************************************************************** Assign geometric reciprocal-temperature schedule **********************************************************************************************/ NUM = 0; TEST_STAT = 0; for (I=1; I<=KR; I++) { RMU[I] = 0; for (J=1; J<=KC; J++) RMU[I] = RMU[I] + X[I][J]; NUM = NUM + RMU[I]; } FNUM = NUM; for (J=1; J<=KC; J++) { CMU[J] = 0; for (I=1; I<=KR; I++) CMU[J] = CMU[J] + X[I][J]; } DUM1 = 0.0; DUM2 = NUM; for(J=1; J<=KC; J++) { DUM1 = DUM1 + sq(CMU[J]); DUM2 = min(DUM2,CMU[J]); } HMAX = sq((NUM - DUM2)) + DUM1 - sq(DUM2); if(OPTIONC == YES) { GMAX = NUM*(NUM/DUM2 - 1); GGMAX = ct_g_max(KR,KC,&RMU[0],&CMU[0]); } for (I=1; I<=KR; I++) { for (J=1; J<=KC; J++) { DUM = RMU[I]*CMU[J]/FNUM; TEST_STAT = TEST_STAT + sq((X[I][J] - DUM))/DUM; } } /******************read initial table from disk******************************/ if (READ == YES) for (I=1; I<=KR; I++) for (J=1; J<=KC; J++) fscanf(fp1, "%15.0f", X[I][J]); /*************************************************************************/ B[1] = 0.0; B[2] = log(1.0 + PSI)/HMAX; BR = -log(1.0 - exp(-log(1.0 + PSI)/KC)); R = 2 + ceil(log(BR/B[2])/log(1.0 + THETA)); B[R] = BR; for (I=3; I<= R-1; I++) B[I] = (1.0+THETA)*B[I-1]; if(OPTIONC == YES && KSTAR > 0 && KSTAR < R) { for(I =0; I<=R-KSTAR-1; I++) B[R-I+1] = B[R-I]; B[KSTAR+1] = 1*B[KSTAR]; R = R +1; } NSTAGES = R; CUM = 0.0; CVAR = 0.0; CRATIO = 0.0; OLDC = 0; for(I=1; I<=NSTAGES; I++) { ZF[I] = 0.0; HMR[I] = 0.0; SUMXI[I] = 0.0; CMIN[I] = 0.0; SUMXIA[I] = 0.0; } /************************************************************************************* Begin independent replication *************************************************************************************/ for(LL=1; LL<=NREPS; LL++) { HOLD = 0.0; SLOG = 0.0; COLD = 0.0; CHLOG = 0.0; CS[LL] = 0.0; if (LL==1) { fprintf(fp, " %d by %d contingency table (Hastings-Metropolis, quadratic penalty and randomly sampling within rows)\n\n", KR,KC); fprintf(fp, " a = %d b = %d m = %10.0f r = %d OPTIONC = %d switch scaling on stage = %d switch nominating kernels on stage = %d\n\n warm up interval= %d t = %15.0f n= %d psi = %6.3f theta = %12.6e\n\n", KR,KC,NUM,R,OPTIONC,KSTAR,SW,WU,T,NREPS,PSI,THETA); } for (J=1; J<=KC; J++) { SC[J] = -CMU[J]; for(I=1; I<=KR; I++) SC[J] = SC[J] + X[I][J]; } for (I=1; I<=KR; I++) { SR[I] = - RMU[I]; KZ[I] = 0; for (J=1; J<=KC; J++) { SR[I] = SR[I] + X[I][J]; if(X[I][J] == 0) KZ[I] = KZ[I] + 1; DUM = RMU[I]*CMU[J]/FNUM; COLD = COLD + sq((X[I][J] - DUM))/DUM; } } /*********************************************************************** Sets up table for sampling row I ***********************************************************************/ QQQ = 0; DELTA = 0; for (I=1; I<=KR; I++) { QQQ = QQQ + KC - KZ[I]; /* QQQ = # of nonzero entries in table */ for (J=1; J<=KC-KZ[I]; J++) W[DELTA + J] = I; DELTA = DELTA + KC - KZ[I]; } for (I=1; I<=KR; I++) for(J=1; J<=KC; J++) Y[I][J] = X[I][J]; /************************************************************************ Quadratic penalty function ************************************************************************/ for(J=1; J<=KC; J++) HOLD = HOLD + SC[J]*SC[J]; /************************************************************************ Begin stage **************************************************************************/ for (STAGE = NSTAGES; STAGE>=LSTAGE; STAGE--) /* (i+1,i) strategy */ { SUMB = 0.0; SUM = 0.0; SSL = 0.0; SUC = 0.0; SUD = 0.0; /************************************************************************ Option C_k scaling *************************************************************************/ DUM3 = (NSTAGES/STAGE - 1)/(NSTAGES - 1); GAM = GGMAX*DUM3 + TEST_STAT*(1 - DUM3); if(STAGE > LSTAGE) { DUM4 = (NSTAGES/(STAGE -1) - 1)/(NSTAGES -1); GAM2 = GGMAX*DUM4 + TEST_STAT*(1 - DUM4); } if(STAGE >= KSTAR) GAM = TEST_STAT; if(STAGE > KSTAR) GAM2 = TEST_STAT; /************************************************************************* Begin sample path *************************************************************************/ for(STEP=1; STEP<=T+WU; STEP++) { CCC: U = genrand_real3(); PICK = ceil(U*QQQ); RI = W[PICK]; IP = RI; AAA: U = genrand_real3(); JP = ceil(KC*U); if (Y[IP][JP] == 0) goto AAA; U = genrand_real3(); RJ = ceil(U*(KC - 1)); if (RJ >= JP) RJ = RJ + 1; PPP = Y[RI][RJ] + Y[IP][JP]; DUM1 = RMU[RI]*CMU[RJ]/NUM; DUM2 = RMU[IP]*CMU[JP]/NUM; INCX = 1; INCY = 1; if(OPTIONC == NO || STAGE == LSTAGE) { if (STAGE >= SW) { ZZ = Y[RI][RJ] + 1; } else { U = genrand_real3(); ZZ = floor(U*PPP); if (ZZ >= Y[RI][RJ]) ZZ = ZZ + 1; } } else { AX = PPP*DUM1/(DUM1 + DUM2) - Y[RI][RJ]; BX = DUM1*DUM2/(DUM1 + DUM2)*(GAM2 - COLD); RADX = sqrt(sq((AX)) + BX); LX = max(0,ceil(Y[RI][RJ] + AX - RADX)); UX = min(PPP,floor(Y[RI][RJ] + AX + RADX)); INCX = UX - LX + 1; if((Y[RI][RJ] - LX)*(UX - Y[RI][RJ]) >= 0) INCX = INCX - 1; U = genrand_real3(); ZZ = LX + floor(U*INCX); if(ZZ >= Y[RI][RJ]) ZZ = ZZ + 1; } CNEW = COLD + (sq((ZZ - DUM1)) - sq((Y[RI][RJ] - DUM1)))/DUM1 + (sq((PPP - ZZ -DUM2)) - sq((Y[IP][JP] - DUM2)))/DUM2; if(OPTIONC == YES && STAGE > LSTAGE) { AY = PPP*DUM1/(DUM1 + DUM2) - ZZ; BY = DUM1*DUM2/(DUM1 + DUM2)*(GAM2 - CNEW); RADY = sqrt(sq((AY)) + BY); LY = max(0,ceil(ZZ + AY - RADY)); UY = min(PPP,floor(ZZ + AY + RADY)); INCY = UY - LY + 1; if((ZZ -LY)*(UY - ZZ) >= 0) INCY = INCY -1; } BIGQ = 1.0; if (STAGE < SW) { if (Y[RI][RJ] > 0) BIGQ = BIGQ/2.0; if (ZZ > 0 && ZZ < PPP) BIGQ = 2.0*BIGQ; } RR = 0; /* compute RR=# of zeros that become non zeros */ if (Y[RI][RJ] == 0) RR = RR + 1; if (ZZ == 0) RR = RR - 1; if (ZZ == PPP) RR = RR - 1; BIGQ = BIGQ*QQQ/(QQQ + RR); BIGQ = BIGQ*INCX/INCY; HNEW = HOLD; /**********************Quadratic penalty function update******************************/ HNEW = HNEW + 2.0*(ZZ - Y[RI][RJ])*(ZZ - Y[RI][RJ] + SC[RJ] - SC[JP]); /***************************************************************************************/ U = genrand_real3(); if(U <= BIGQ*exp(-B[STAGE]*(HNEW-HOLD))) { HOLD = HNEW; if(STEP > WU ) HMR[STAGE] = HMR[STAGE] + 1.0; if (RR > 0) { KZ[RI] = KZ[RI] - 1; QQQ = QQQ + 1; W[QQQ] = RI; } if (RR < 0) { KZ[RI] = KZ[RI] + 1; W[PICK]= W[QQQ]; QQQ = QQQ - 1; } /*********************************************************************** Update ************************************************************************/ DUM1 = RMU[RI]*CMU[RJ]/NUM; DUM2 = RMU[IP]*CMU[JP]/NUM; COLD = COLD + (sq((ZZ - DUM1)) - sq((Y[RI][RJ] - DUM1)))/DUM1 + (sq((PPP - ZZ -DUM2)) - sq((Y[IP][JP] - DUM2)))/DUM2; SC[RJ] = SC[RJ] - Y[RI][RJ] + ZZ ; SC[JP] = SC[JP] - Y[IP][JP] + PPP - ZZ; Y[RI][RJ] = ZZ; Y[IP][JP] = PPP - ZZ; if (HOLD == 0.0) { if(STAGE > LSTAGE) { if(OPTIONC == NO) { for (I=1; I<=KR; I++) for (J=1; J<=KC; J++) X[I][J] = Y[I][J]; } else { if (COLD <= TEST_STAT) for (I=1; I<=KR; I++) for (J=1; J<=KC; J++) X[I][J] = Y[I][J]; } } } } /************************************************************************ If warm-up interval is over, begin data collection ************************************************************************/ if(STEP > WU) { OLDC =max(OLDC,COLD); if (HOLD == 0.0) CMIN[STAGE] = CMIN[STAGE] +1.0; if(STAGE == NSTAGES) { DDD = 1.0; if(HOLD > 0.0) DDD = 0.0; if(COLD <= GAM) { CUM = CUM + DDD; CS[LL] = CS[LL] + DDD; } } else { DDD = exp(-(B[STAGE+1]-B[STAGE])*HOLD); } if (OPTIONC == YES && COLD > GAM) DDD = 0.0; if (COLD <= GAM) SUMXIA[STAGE] += sq(DDD); SUMB = SUMB + DDD; SUMXI[STAGE] = SUMXI[STAGE] + sq(DDD); if(COLD <= GAM) SUD = SUD + DDD; if (STAGE == LSTAGE) { SUC = SUC + 1.0; } else { if(COLD <= GAM2) SUC = SUC + 1.0; } } } /******************************************************************************* Begin data collection on this independent replication ********************************************************************************/ A[STAGE][LL] = SUMB; STN[STAGE][LL] = SUD; STD[STAGE][LL] = SUC; } } /******************************************************************************* end STAGE ************************ end independent replication***************************/ if(OPTIONC == NO) { fprintf(fp, " Option A\n"); fprintf(fp, " I B[I] LAMBDA s.e.(LAMBDA) r.e.(LAMBDA) XI RHO HM RATE r.e.(GAMMA) U ALPHA\n *************************************\n"); } if(OPTIONC == YES) { fprintf(fp, " Option C\n*******************************************************************************************\n"); fprintf(fp, " I B[I] GAMMA' s.e.(GAMMA') r.e.(GAMMA') U' ALPHA' HM rate\n"); } DOT = 0.0; C2DOT = 0.0; E2DOT = 0.0; for(I=LSTAGE; I<=NSTAGES; I++) { SA = 0.0; SG = NREPS*T; SN = 0.0; SD = 0.0; for(J=1; J<=NREPS; J++) { SA += A[I][J]; SN += STN[I][J]; if(I == LSTAGE) STD[I][J] = T; SD += STD[I][J]; } LAM = SA/SG; LAMA = SN/SD; SLOG += log(LAM); CHLOG +=log(SN/SD); XI = SUMXI[I]/LAM/LAM/(SG - 1.0) - SG/(SG - 1.0); XIA = (SUMXIA[I]/SN/SN - 1/SD)*SG; HMR[I] = HMR[I]/SG; ZF[I] = ZF[I]/SG; S2A = 0.0; DS2A = 0.0; ES2A = 0.0; for(J=1; J<=NREPS; J++) { DA = A[I][J]/T - LAM; S2A = S2A + sq(DA); if(I > LSTAGE) { DE = STN[I][J]/(SN/NREPS) - STD[I][J]/(SD/NREPS); EE = DE - A[I][J]/T/LAM + 1.0; } else { DE = STN[I][J]/(SN/NREPS) - 1; EE = DE; } DS2A = DS2A + sq(DE); ES2A = ES2A + sq(EE); } S2A = S2A/(NREPS - 1.0); VAR = S2A/NREPS; RHO[I] = VAR*SG/LAM/LAM/XI; REL = 0.0; if (LAM > 0.0) REL = sqrt(VAR)/LAM; DS2A = DS2A/(NREPS - 1.0)/NREPS; C2REL = sqrt(DS2A); C2DOT = C2DOT + log(1.0 + sq(C2REL)); ES2A = ES2A/(NREPS - 1.0)/NREPS; E2REL = sqrt(ES2A); E2DOT += log(1.0 + sq(E2REL)); if(I < 10) fprintf(fp, " "); if(OPTIONC == NO) fprintf(fp, "%i %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %9.2e *%12.3e %12.3e %12.3e\n", I,B[I],LAM,LAM*REL,REL,XI,RHO[I],HMR[I],C2REL,XIA,SG*DS2A/XIA); if(OPTIONC == YES) fprintf(fp, "%i %12.3e %12.3e %12.3e %12.3e %12.3e %12.3e %9.2e\n", I,B[I],LAM,C2REL,XIA,SG*DS2A/XIA,RHO[I],HMR[I]); DOT = DOT + log(1.0 + REL*REL); } DUM1 = 0.0; AZ = 0; BZ = 0; for (I=1; I<=KR; I++) { DUM1 = log(KC); DUM2 = log(1.0); Z1 = exp(DUM1 + SLOG/KR); Z2 = exp(DUM1 + CHLOG/KR); MM = min(KC,RMU[I]); for (J=2; J<=MM; J++) { DUM1 = DUM1 + log((KC - J + 1.0)/J); DUM2 = DUM2 + log((RMU[I] - J + 1.0)/(J -1.0)); Z1 = Z1 + exp(DUM1 + DUM2 + SLOG/KR); Z2 = Z2 + exp(DUM1 + DUM2 + CHLOG/KR); } AZ = AZ + log(Z1); BZ = BZ + log(Z2); } /******************************************************************************************************** total count ********************************************************************************************************/ Z = exp(AZ); CHZ = exp(BZ); REZ = sqrt(exp(DOT) - 1.0); SEZ = Z*REZ; /********************************************************************************************************** Option A: chi squared analysis **********************************************************************************************************/ C2REZ = sqrt(exp(C2DOT) -1.0); /* count < test statistic */ C2SEZ = CHZ*C2REZ; E2REZ = sqrt(exp(E2DOT) -1.0); /* proportion */ E2SEZ = (CHZ/Z)*E2REZ; /******************************************************************************************************** Option B: chi squared statistic analysis ********************************************************************************************************/ I = NSTAGES; for(J=1; J<=NREPS; J++) { DS = CS[J]/T - CUM/SG; CVAR += sq(DS); DA = A[I][J]/T - LAM; ND = DS/(CUM/SG) - DA/LAM; CRATIO += ND*ND; } CVAR = CVAR/(NREPS - 1.0)/NREPS; CRATIO = CRATIO/(NREPS - 1.0)/NREPS; CREL = sqrt(CVAR)/(CUM/SG); CHIZ = Z*CUM/SG/LAM; /* count < test statistic */ DOT = DOT - log(1.0 + sq(REL)) + log(1.0 + sq(CREL)); CREZ = sqrt(exp(DOT) - 1.0); CSEZ = CHIZ*CREZ; AC = CUM/SG/LAM; /* proportion */ /******************************************************************************************************** printed output ********************************************************************************************************/ fprintf(fp, "\n log_10(|X|)= %15.6e max H = %15.6e\n\n", log10(Z)-log10(exp(SLOG)),HMAX); if(OPTIONC == YES) fprintf(fp, "\n ~g_max(A) = %15.6e max observed g(x) = %15.6e\n\n", GGMAX,OLDC); fprintf(fp, "*******************************************************************************************************************\n"); fprintf(fp, "*******************************************************************************************************************\n"); if(OPTIONC == NO) { fprintf(fp, "NUMBER OF CONTINGENCY TABLES\n\n"); fprintf(fp, " |A|= %12.6e S.E.= %12.6e R.E. = %12.6e\n", Z,SEZ,REZ); fprintf(fp, "*******************************************************************************************************************\n"); fprintf(fp, "NUMBER WITH CHI-SQUARED STATISTIC <= g(x*) =%12.6e\n\n", TEST_STAT); fprintf(fp, "*******************************************************************************************************************\n"); fprintf(fp, "Option A: |A(C)| = %12.6e S.E. = %12.6e R.E. = %12.6e\n\n", CHZ, C2SEZ,C2REZ); fprintf(fp, " proportion: C := [0,%12.6e) |A(C)|/|A| = %12.6e S.E. = %12.6e R.E. = %12.6e\n", TEST_STAT,CHZ/Z,E2SEZ,E2REZ); fprintf(fp, "*******************************************************************************************************************\n"); fprintf(fp, "Option B: |A(C)| = %12.6e S.E. = %12.6e R.E. = %12.6e\n\n", CHIZ,CSEZ,CREZ); fprintf(fp, " proportion: C := [0,%12.6e) |A(C)|/|A| = %12.6e S.E. = %12.6e R.E. = %12.6e\n", TEST_STAT,AC,AC*sqrt(CRATIO),sqrt(CRATIO)); fprintf(fp, "*******************************************************************************************************************\n"); } else { fprintf(fp, "Option C: NUMBER OF CONTINGENCY TABLES <= %15.6e\n\n", TEST_STAT); fprintf(fp, " |A|= %12.6e S.E.= %12.6e R.E. = %12.6e\n", Z,SEZ,REZ); } fprintf(fp, "*******************************************************************************************************************\n\n"); /****************************print pr[X in A]***************************************/ fprintf(fp, " STAGE pr[X in A]\n"); for(I=1; I<=NSTAGES; I++) { if(I < 10)fprintf(fp, " "); fprintf(fp, "%d %12.4e\n", I,CMIN[I]/NREPS/T); } /*******************Save last contingency table for next initialization***************/ fclose(fp1); fp1 = fopen("ct.hm.initial.txt", "w"); for (I=1; I<=KR; I++) for(J=1; J<=KC; J++) fprintf(fp1, "%10.0f", X[I][J]); /*****************save last random numbers as seed for next initialization***********/ fsave_seeds(); /************************************************************************************/ fclose(fp); fclose(fp1); } /* end ct.chisq.c procedure */ double ct_g_max(int KR, int KC, double MR[], double MC[]) { /* begin ct_g_max.c procedure */ int I,II,J,K,L,W; double X[20][20],Y[20][20]; double Z[20][20],RSUM[20],CSUM[20]; double CHISQ,NEWCHISQ=0; double MM; double Q; double SUM,SUMA,SUMB; double A=0,AMIN,AMAX,B,DUM,DUM1,DUM2; /************************************************************************************************************************** ct_g_max.c orders row sums, then column sums and then computes max_{x in A} g(x), where g(x) is the chi-squared statistic for x. It is called by ct.chisq.c. **************************************************************************************************************************/ MM = 0; for(I=1; I<=KR; I++) MM = MM + MR[I]; for(I=1; I<=KR; I++) RSUM[I] = MR[I]; for(J=1; J<=KC; J++) CSUM[J] = MC[J]; for(I=1; I<=KR; I++) { for(J=1; J<=KC; J++) X[I][J] =0; } SUMA = MM; SUMB = MM; AAA: for(I=1; I<=KR; I++) { for(J=1; J<=KC; J++) { Z[I][J] = 0; if(RSUM[I]*CSUM[J] > 0) Z[I][J] = min((RSUM[I]/CSUM[J]),(CSUM[J]/RSUM[I])); } } Q = 0; for(I=1; I<=KR; I++) { for(J=1; J<=KC; J++) { if (Q < Z[I][J]) { Q = Z[I][J]; K = I; L = J; } } } if(Q == RSUM[K]/CSUM[L]) { X[K][L] = RSUM[K]; CSUM[L] = CSUM[L] - RSUM[K]; RSUM[K] = 0; } else { X[K][L] = CSUM[L]; RSUM[K] = RSUM[K] - CSUM[L]; CSUM[L] = 0; } SUMA = SUMA - X[K][L]; SUMB = SUMB - X[K][L]; SUM = SUM + X[K][L]; if( SUM < MM ) goto AAA; CHISQ = 0; for(I=1; I<=KR; I++) { for(J=1; J<=KC; J++) CHISQ = CHISQ + X[I][J]*X[I][J]/MR[I]/MC[J]; } CHISQ = MM*(CHISQ -1); /************************************************************************************************* Incremental improvemnts *************************************************************************************************/ BBB: II = 0; for(I=1; I<=KR-1; I++) { for(J=1; J<=KC-1; J++) { for(K=I+1; K<=KR; K++) { for(L=J+1; L<=KC; L++) { B = -2*(X[I][J]/MR[I]/MC[J] - X[I][L]/MR[I]/MC[L] -X[K][J]/MR[K]/MC[J] + X[K][L]/MR[K]/MC[L])/ (1/MR[I]/MC[J] + 1/MR[I]/MC[L] + 1/MR[K]/MC[J] + 1/MR[K]/MC[L]); AMIN = -min(X[I][J],X[K][L]); AMAX = min(X[I][L],X[K][J]); if (B < 0) { if (AMIN < B) A = AMIN; } if(B > 0) { if (AMAX > B) A = AMAX; } /*************************************************************************************************** For B= 0, A=0 which says no change or as follows Question: In this case are AMIN and AMAX =0? */ if(B == 0) { A = AMIN; if(AMAX > abs(AMIN)) A = AMAX; } /***************************************************************************************************/ if (A > 0 || A < 0) { II = 1; CHISQ = CHISQ + MM*A*((2*X[I][J] + A)/MR[I]/MC[J] - (2*X[I][L] - A)/MR[I]/MC[L] - (2*X[K][J] - A)/MR[K]/MC[J] + (2*X[K][L] + A)/MR[K]/MC[L]); X[I][J] = X[I][J] + A; X[I][L] = X[I][L] - A; X[K][J] = X[K][J] - A; X[K][L] = X[K][L] + A; A = 0; } } } } } if (II >0 ) goto BBB; return CHISQ; } /* end ct_g_max.c procedure */