/* The following C program finds the value of A (specifically for our frog example) used in eigenvalue comparison in Diaconis and Saloff-Coste (1993), "Comparison theorems for reversible Markov chains." Ann. Appl. Prob. 3, 696-730 */ #define NODES 200 /* the size of the state space */ #define M1 10 /* P1 has edges between nodes within M1 of each other */ #define M2 9 /* P2 has edges between nodes within M2 of each other */ #define BETA 10 /* the target distribution is proportional to (dist(x,0)+1)^-BETA where dist(x,0) is the distance from node x to node 0 */ #include #include #include int findLength(int, int, int, int); int distance (int); int main() { int min, i, j, z, w, x, y, length; int searchstart, searchstart2, searchend, searchend2; double alpha, max=0, sum, denom=0; double P1[NODES][NODES], P2[NODES][NODES], pi1[NODES], pi2[NODES]; for (i=0; i= pi1[i]) alpha = 1; /* alpha is the probability that the proposed state, j, is accepted given that we are in state i; this is the same for P1 and P2 since they have the same stationary distribution. However, the proposals are different. */ else alpha = pi1[j]/pi1[i]; if ((0 < abs(i-j)) && (abs(i-j) < M1+1)) { P1[i][j] = alpha/(2*M1+1); P1[i][i] = P1[i][i] + (1-alpha)/(2*M1+1); } else if (abs(i-j) >= NODES-M1) /* the highest nodes are connected to the lowest nodes */ { P1[i][j] = alpha/(2*M1+1); P1[i][i] = P1[i][i] + (1-alpha)/(2*M1+1); } else if ((M1+1 <= abs(i-j)) && (abs(i-j) < NODES-M1)) P1[i][i] = 0; if ((0 < abs(i-j)) && (abs(i-j) < M2+1)) { P2[i][j] = alpha/(2*M2+1); P2[i][i] = P2[i][i] + (1-alpha)/(2*M2+1); } else if (abs(i-j) >= NODES-M2) { P2[i][j] = alpha/(2*M2+1); P2[i][i] = P2[i][i] + (1-alpha)/(2*M2+1); } else if ((M2+1 <= abs(i-j)) && (abs(i-j) < NODES-M2)) P2[i][j] = 0; } } for (z=0; z searchend) searchend += NODES; /* allows wrapping around from NODES-1 to 0 */ for (i=searchstart; i<=searchend; i++) { w = i % NODES; sum = 0; for (x=0; x searchend2) searchend2 += NODES; for (j=searchstart2; j<=searchend2; j++) { y = j % NODES; length = findLength(z, w, x, y); sum = sum + length*pi2[x]*P2[x][y]; } } sum = sum/(pi1[z]*P1[z][w]); if (sum > max) max = sum; } } printf ("A is %g\n", max); exit(0); } /* returns the path length from x to y if (z,w) was one of the steps on the path, otherwise returns 0 */ int findLength(int z, int w, int x, int y) { int left, right, indicator, rem, min, len; indicator = 0; /* changes to 1 when (z,w) is found to be a step on the path from x to y */ /* the edge comes by moving clockwise, e.g. 1->2*/ if (((x= NODES/2))) { left = x; right = (left + M1) % NODES; while ((abs(left-y) > M1) && (abs(left-y)1 */ else if (((x= NODES/2)) || ((yM1) && (abs(right-y)