## Monday, March 8, 2010

### 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    Cumulativecorrect    Frequency     Percent     Frequency      Percent------------------------------------------------------------      0     4996889       49.97       4996889        49.97      1     5003111       50.03      10000000       100.00Proportion                0.5003ASE                       0.000295% Lower Conf Limit      0.500095% Upper Conf Limit      0.5006Exact Conf Limits95% Lower Conf Limit      0.500095% Upper Conf Limit      0.5006ASE under H0              0.0002Z                         1.9676One-sided Pr >  Z         0.0246Two-sided Pr > |Z|        0.0491Sample 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 + 1000real2 = runif(n) * 1000 + 1000pickenv = (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 testdata:  sum(larger1(1e+07)) and 1e+07 number of successes = 5003996, number of trials = 1e+07, p-value =0.01150alternative 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 `

#### 1 comment:

Ken Kleinman said...

Folks-- We revisit this problem in the next post: http://sas-and-r.blogspot.com/2010/03/example-727-probability-question.html.