Example 7.26: probability question

Here's a surprising problem, from the xkcd blog.

Suppose I choose two (different) real numbers, by any process I choose. Then I select one at random (p= .5) to show Nick. Nick must guess whether the other is smaller or larger. Being right 50% of the time is easy. Can he do better?

Of course, it wouldn't be an interesting question if the answer was no. Randall Munroe (author of the xkcd webcomic) offers a solution as follows:

1. Call the number you see x.  
      Calculate p(x) = (e^x)/(1 + e^x)
2. Draw u, a random uniform(0,1).
3. If u < p(x) then guess the hidden number is 
      smaller than the revealed one.  Otherwise,
      guess that it is smaller.

Munroe shows that the probability of guessing correctly is

N = 0.5 * (1-p(A)) + 0.5 * p(B)

where A is the smaller of the two numbers. In practice, this only gains us an advantage when the observed value is between about -20 and 20-- outside that range, calculation of p(x) results in numbers so close to 0 or 1 that no improvement over a 50% rate can be achieved. If I knew Nick were using this method, I could keep him to a 50% rate by always choosing numbers over 1000 or so. He can improve Munroe's solution by using p(x/100000), though at the cost of a poorer improvement near 0.

We'll check the logic by simulating the experiment.

SAS

The SAS solution requires a lot of looping (section 1.11.1) and conditional execution (section 1.11.2).

data test;
  do i = 1 to 10000000;
    real1 = uniform(0) * 1000 + 1000;
    real2 = uniform(0) * 1000 + 1000;
 if uniform(0) gt .5 then do;
      env1 = real1;
   env2 = real2;
   end;
 else do;
   env1 = real2;
   env2 = real1;
   end;
 if uniform(0) lt (1/(1 + exp(-1 * env1/100000))) then guess = "l";
 else guess = "h";
 correct = ((env1 < env2) and (guess eq "h")) or 
              ((env1 > env2) and (guess eq "l")) ;
 output;
 end;
run;

We can test the performance using proc freq with the binomial option (section 2.1.9).

proc freq data = test;
tables correct / binomial;
run;

The FREQ Procedure

                                    Cumulative    Cumulative
correct    Frequency     Percent     Frequency      Percent
------------------------------------------------------------
      0     4996889       49.97       4996889        49.97
      1     5003111       50.03      10000000       100.00


Proportion                0.5003
ASE                       0.0002
95% Lower Conf Limit      0.5000
95% Upper Conf Limit      0.5006

Exact Conf Limits
95% Lower Conf Limit      0.5000
95% Upper Conf Limit      0.5006

ASE under H0              0.0002
Z                         1.9676
One-sided Pr >  Z         0.0246
Two-sided Pr > |Z|        0.0491

Sample Size = 10000000

We can just barely detect the improvement-- but it took 10,000,000 trials to find CI that reliably exclude the possibility of a 50% success rate.




R

We begin by writing a function to calculate p(x/100000).

xkcdprob = function(x) {
1/(1 + exp(-1 * x/100000))}

Then we write a function to perform the experiment n times. Here we use the same process for picking the two numbers, choosing them to be Uniform between 1000 and 2000. In the code below we use the ifelse() function (section 1.11.2) to decide which envelope is shown first. The vectorization in R allows us to avoid loops entirely in writing the code.

larger1 = function(n) {
real1 = runif(n) * 1000 + 1000
real2 = runif(n) * 1000 + 1000
pickenv = (runif(n) < .5)
env1 = ifelse(pickenv,real1,real2)
env2 = ifelse(!pickenv,real1,real2)
guess = ifelse(runif(n) < xkcdprob(env1),"lower","higher")
correct =  ((env1 < env2) & (guess == "higher"))  | ((env1 > env2) &
   (guess == "lower"))
return(correct)
}

Then we can run the experiment and check its results in one nested call, using the binom.test() function (section 2.1.9) to see whether the observed proportion is different from 0.5.

> binom.test(sum(larger1(10000000)),10000000)

        Exact binomial test

data:  sum(larger1(1e+07)) and 1e+07 
number of successes = 5003996, number of trials = 1e+07, p-value =
0.01150
alternative hypothesis: true probability of success is not equal to 0.5 
95 percent confidence interval:
 0.5000897 0.5007095 
sample estimates:
probability of success 
             0.5003996