406 probability/flips/waiting.time.p
Description
This article is from the Puzzles FAQ, by Chris Cole chris@questrel.questrel.com and Matthew Daly mwdaly@pobox.com with numerous contributions by others.
406 probability/flips/waiting.time.p
Compute the expected waiting time for a sequence of coin flips, or the
probabilty that one sequence of coin flips will occur before another.
probability/flips/waiting.time.s
Here's a C program I had lying around that is relevant to the
current discussion of coin-flipping. The algorithm is N^3 (for N flips)
but it could certainly be improved. Compile with
cc -o flip flip.c -lm
-- Guy
_________________ Cut here ___________________
#include <stdio.h>
#include <math.h>
char *progname; /* Program name */
#define NOT(c) (('H' + 'T') - (c))
/* flip.c -- a program to compute the expected waiting time for a sequence
of coin flips, or the probabilty that one sequence
of coin flips will occur before another.
Guy Jacobson, 11/1/90
*/
main (ac, av) int ac; char **av;
{
char *f1, *f2, *parseflips ();
double compute ();
progname = av[0];
if (ac == 2) {
f1 = parseflips (av[1]);
printf ("Expected number of flips until %s = %.1f\n",
f1, compute (f1, NULL));
}
else if (ac == 3) {
f1 = parseflips (av[1]);
f2 = parseflips (av[2]);
if (strcmp (f1, f2) == 0) {
printf ("Can't use the same flip sequence.\n");
exit (1);
}
printf ("Probability of flipping %s before %s = %.1f%%\n",
av[1], av[2], compute (f1, f2) * 100.0);
}
else
usage ();
}
char *parseflips (s) char *s;
{
char *f = s;
while (*s)
if (*s == 'H' || *s == 'h')
*s++ = 'H';
else if (*s == 'T' || *s == 't')
*s++ = 'T';
else
usage ();
return f;
}
usage ()
{
printf ("usage: %s {HT}^n\n", progname);
printf ("\tto get the expected waiting time, or\n");
printf ("usage: %s s1 s2\n\t(where s1, s2 in {HT}^n for some fixed n)\n",
progname);
printf ("\tto get the probability that s1 will occur before s2\n");
exit (1);
}
/*
compute -- if f2 is non-null, compute the probability that flip
sequence f1 will occur before f2. With null f2, compute
the expected waiting time until f1 is flipped
technique:
Build a DFA to recognize (H+T)*f1 [or (H+T)*(f1+f2) when f2
is non-null]. Randomly flipping coins is a Markov process on the
graph of this DFA. We can solve for the probability that f1 precedes
f2 or the expected waiting time for f1 by setting up a linear system
of equations relating the values of these unknowns starting from each
state of the DFA. Solve this linear system by Gaussian Elimination.
*/
typedef struct state {
char *s; /* pointer to substring string matched */
int len; /* length of substring matched */
int backup; /* number of one of the two next states */
} state;
double compute (f1, f2) char *f1, *f2;
{
double solvex0 ();
int i, j, n1, n;
state *dfa;
int nstates;
char *malloc ();
n = n1 = strlen (f1);
if (f2)
n += strlen (f2); /* n + 1 states in the DFA */
dfa = (state *) malloc ((unsigned) ((n + 1) * sizeof (state)));
if (!dfa) {
printf ("Ouch, out of memory!\n");
exit (1);
}
/* set up the backbone of the DFA */
for (i = 0; i <= n; i++) {
dfa[i].s = (i <= n1) ? f1 : f2;
dfa[i].len = (i <= n1) ? i : i - n1;
}
/* for i not a final state, one next state of i is simply
i+1 (this corresponds to another matching character of dfs[i].s
The other next state (the backup state) is now computed.
It is the state whose substring matches the longest suffix
with the last character changed */
for (i = 0; i <= n; i++) {
dfa[i].backup = 0;
for (j = 1; j <= n; j++)
if ((dfa[j].len > dfa[dfa[i].backup].len)
&& dfa[i].s[dfa[i].len] == NOT (dfa[j].s[dfa[j].len - 1])
&& strncmp (dfa[j].s, dfa[i].s + dfa[i].len - dfa[j].len + 1,
dfa[j].len - 1) == 0)
dfa[i].backup = j;
}
/* our dfa has n + 1 states, so build a system n + 1 equations
in n + 1 unknowns */
eqsystem (n + 1);
for (i = 0; i < n; i++)
if (i == n1)
equation (1.0, n1, 0.0, 0, 0.0, 0, -1.0);
else
equation (1.0, i, -0.5, i + 1, -0.5, dfa[i].backup, f2 ? 0.0 : -1.0);
equation (1.0, n, 0.0, 0, 0.0, 0, 0.0);
free (dfa);
return solvex0 ();
}
/* a simple gaussian elimination equation solver */
double *m, **M;
int rank;
int neq = 0;
/* create an n by n system of linear equations. allocate space
for the matrix m, filled with zeroes and the dope vector M */
eqsystem (n) int n;
{
char *calloc ();
int i;
m = (double *) calloc (n * (n + 1), sizeof (double));
M = (double **) calloc (n, sizeof (double *));
if (!m || !M) {
printf ("Ouch, out of memory!\n");
exit (1);
}
for (i = 0; i < n; i++)
M[i] = &m[i * (n + 1)];
rank = n;
neq = 0;
}
/* add a new equation a * x_na + b * x_nb + c * x_nc + d = 0.0
(note that na, nb, and nc are not necessarily all distinct.) */
equation (a, na, b, nb, c, nc, d) double a, b, c, d; int na, nb, nc;
{
double *eq = M[neq++]; /* each row is an equation */
eq[na + 1] += a;
eq[nb + 1] += b;
eq[nc + 1] += c;
eq[0] = d; /* column zero holds the constant term */
}
/* solve for the value of variable x_0. This will go nuts if
therer are errors (for example, if m is singular) */
double solvex0 ()
{
register i, j, jmax, k;
register double max, val;
register double *maxrow, *row;
for (i = rank; i > 0; --i) { /* for each variable */
/* find pivot element--largest value in ith column*/
max = 0.0;
for (j = 0; j < i; j++)
if (fabs (M[j][i]) > fabs (max)) {
max = M[j][i];
jmax = j;
}
/* swap pivot row with last row using dope vectors */
maxrow = M[jmax];
M[jmax] = M[i - 1];
M[i - 1] = maxrow;
/* normalize pivot row */
max = 1.0 / max;
for (k = 0; k <= i; k++)
maxrow[k] *= max;
/* now eliminate variable i by subtracting multiples of pivot row */
for (j = 0; j < i - 1; j++) {
row = M[j];
if (val = row[i]) /* if variable i is in this eq */
for (k = 0; k <= i; k++)
row[k] -= maxrow[k] * val;
}
}
/* the value of x0 is now in constant column of first row
we only need x0, so no need to back-substitute */
val = -M[0][0];
free (M);
free (m);
return val;
}
_________________________________________________________________
Guy Jacobson (201) 582-6558 AT&T Bell Laboratories
uucp: {att,ucbvax}!ulysses!guy 600 Mountain Avenue
internet: guy@ulysses.att.com Murray Hill NJ, 07974
Continue to:
My Books
