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

## 1 comment:

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

Post a Comment