Example 8.35: Grab true (not pseudo) random numbers; passing API URLs to functions or macros

Usually, we're content to use a pseudo-random number generator. But sometimes we may want numbers that are actually random-- an example might be for randomizing treatment status in a randomized controlled trial.

The site Random.org provides truly random numbers based on radio static. For long simulations, its quota system may prevent its use. But for small to moderate needs, it can be used to provide truly random numbers. In addition, you can purchase larger quotas if need be.

The site provides APIs for several types of information. We'll write functions to use these to pull vectors of uniform (0,1) random numbers (of 10^(-9) precision) and to check the quota. To generate random variates from other distributions, you can use the inverse probability integral transform (section 1.10.8).

The coding challenge here comes in integrating quotation marks and special characters with function and macro calls.

SAS
In SAS, the challenging bit is to pass the desired number of random numbers off to the API, though the macro system. This is hard because the API includes the special characters ?, ", and especially &. The ampersand is used by the macro system to denote the start of a macro variable, and is used in APIs to indicate that an additional parameter follows.

To avoid processing these characters as part of the macro syntax, we have to enclose them within the macro quoting function %nrstr. We use this approach twice, for the fixed pieces of the API, and between them insert the macro variable that contains the number of random numbers desired. Also note that the sequence %" is used to produce the quotation mark. Then, to unmask the resulting character string and use it as intended, we %unquote it. Note that the line break shown in the filename statement must be removed for the code to work.

Finally, we read data from the URL (section 1.1.6) and transform the data to take values between 0 and 1.

%macro rands (outds=ds, nrands=);
filename randsite url %unquote(%nrstr(%"http://www.random.org/integers/?num=)
 &nrands%nrstr(&min=0&max=1000000000&col=1&base=10&format=plain&rnd=new%"));
proc import datafile=randsite out = &outds dbms = dlm replace; 
  getnames = no; 
run;

data &outds;
set &outds;
  var1 = var1 / 1000000000;
run;
%mend rands;

/* an example macro call */
%rands(nrands=25, outds=myrs);

The companion macro to find the quota is slightly simpler, since we don't need to insert the number of random numbers in the middle of the URL. Here, we show the quota in the SAS log; the file print syntax, shown in Example 8.34, can be used to send it to the output instead.

%macro quotacheck;
filename randsite url %unquote(%nrstr(%"http://www.random.org/quota/?format=plain%"));
proc import datafile=randsite out = __qc dbms = dlm replace; 
getnames = no; 
run;

data _null_;
set __qc;
put "Remaining quota is " var1 "bytes";
run;
%mend quotacheck;

/* an example macro call */
%quotacheck;

R

Two R functions are shown below. While the problem isn't as difficult as in SAS, it is necessary to enclose the character string for the URL in the as.character() function (section 1.4.1).

truerand = function(numrand) {
  read.table(as.character(paste("http://www.random.org/integers/?num=", 
  numrand, "&min=0&max=1000000000&col=1&base=10&format=plain&rnd=new", 
  sep="")))/1000000000
}

quotacheck = function() {
  line = as.numeric(readLines("http://www.random.org/quota/?format=plain"))
  return(line)
}

Example 7.23: the Monty Hall problem

The Monty Hall problem illustrates a simple setting where intuition often leads to a solution different from formal reasoning. The situation is based on the game show Let's Make a Deal. First, Monty puts a prize behind one of three doors. Then the player chooses a door. Next, (without moving the prize) Monty opens an unselected door, revealing that the prize is not behind it. The player may then switch to the other non-selected door. Should the player switch?

Many people see that there are now two doors to choose between, and feel that since Monty can always open a non-prize door, there's still equal probability for each door. If that were the case, the player might as well keep the original door.

This approach is so attractive that when Marilyn Vos Savant asserted that the player should switch (in her Parade magazine column), there were reportedly 10,000 letters asserting she was wrong. There was also an article published in The American Statistician. As reported on the Wikipedia page, out of 228 subjects in one study, only 13% chose to switch (Granberg and Brown, Personality and Social Psychology Bulletin 21(7): 711-729). In her book, The Power of Logical Thinking, vos Savant quotes cognitive psychologist Massimo Piattelli-Palmarini as saying "... no other statistical puzzle comes so close to fooling all the people all the time" and "that even Nobel physicists systematically give the wrong answer, and that they insist on it, and they are ready to berate in print those who propose the right answer."

One correct intuitive route is to observe that Monty's door is fixed. The probability that the player has the right door is 1/3 before Monty opens the non-prize door, and remains 1/3 after that door is open. This means that the probability the prize is behind one of the other doors is 2/3, both before and after Monty opens the non-prize door. After Monty opens the non-prize door, the player gets a 2/3 chance of winning by switching to the remaining door. If the player wants to win, they should switch doors.

The excellent Wikipedia entry referenced above provides additional intuitive tools, as well as variants and history.

One way to prove to yourself that switch is best is through simulation. In fact, even deciding how to code the problem may be enough to convince yourself to switch.

SAS
In SAS, we assign the prize to a door, then make an initial guess. If the guess was right, Monty can open either door. We'll switch to the other door. Rather than have Monty choose a door, we'll choose one, under the assumption that Monty opened the other one. We do this with a do until loop (section 1.11.1). If our initial guess was wrong, Monty will open the only remaining non-prize door, and when we switch we'll be choosing the prize door.

data mh;
do i = 1 to 100000;
  prize = rand("TABLE",.333,.333);
* Monty puts the prize behind a random door;
  initial_guess = rand("TABLE",.333,.333);
* We make a random initial guess;
* if the initial guess is right, Monte can 
                  open either of the others;
* which means that player can switch to either 
                  of the other doors;
  if initial_guess eq prize then do;
    new_guess = initial_guess;
    * choose a door until it's different from
                  the initial guess-- that's
                  the door we switch to;
    do until (new_guess ne initial_guess);
      new_guess = rand("TABLE",.333,.333);
   end;
    end;
* If the initial guess was wrong, Monte can rule 
          out only one of the other doors;
* which means that we must switch to the right door;
  if initial_guess ne prize then new_guess = prize;
  output;
  end;
run;

To see what happened, we'll summarize the resulting data set. If our initial guess was right, we'll win by keeping it. If switching leads to a win, the new guess is where the prize is. We create new variables to indicate winning using a logical test (as in section 1.4.9).

data mh2;
set mh;
win_by_keep = (initial_guess eq prize);
win_by_switch = (new_guess eq prize);
run;

proc means data = mh2 mean;
var win_by_keep win_by_switch;
run;

The MEANS Procedure

Variable                 Mean
-----------------------------
win_by_keep         0.3332700
win_by_switch       0.6667300
-----------------------------

Note that these two values sum to exactly one. If we chose the prize initially, we always win by keeping it, and if we did not choose the prize initially, we always win by switching. In any event, the simulation supports the conclusion that we should switch.



R

In R, we write two functions. In one, Monty opens a door, choosing at random among the non-chosen doors if the initial choice was correct, or choosing the one non-selected non-prize door if the initial choice was wrong. The other function returns the door chosen by swapping. We use the sample() function (section 1.5.2) to randomly pick one value. We then use these functions on each trial with the apply() statement (section B.5.3).

numsim = 100000
doors = 1:3
opendoor = function(x) { 
 if (x[1]==x[2])
    return(sample(doors[-c(x[1])], 1)) 
 else return(doors[-c(x[1],x[2])])
}
swapdoor = function(x) { return(doors[-c(x[1], x[2])]) }
winner = sample(doors, numsim, replace=TRUE)
choice = sample(doors, numsim, replace=TRUE)
open = apply(cbind(winner, choice), 1, opendoor)
newchoice = apply(cbind(open, choice), 1, swapdoor)

To display the results, we use the cat() statement to show how text and variables can be integrated.

> cat("without switching, won ",
+ round(sum(winner==choice)/numsim*100,1)," 
+ percent of the time.\n", sep="")
without switching, won 33.2 percent of the time.
> cat("with switching, won ",
+ round(sum(winner==newchoice)/numsim*100,1)," 
+ percent of the time.\n", sep="")
with switching, won 66.8 percent of the time.

Example 7.21: Write a function to simulate categorical data

In example 7.20, we showed how to simulate categorical data. But we might anticipate needing to do that frequently. If a SAS function weren't built in and an equivalent R function not available in a package, we could build them from scratch.

SAS

The SAS code is particularly tortured, since we must parse the parameter string to extract the table probabilities. We do this with the count and find (section 1.4.6) and substr (section 1.4.3) functions. We also turn the characters into numbers (section 1.4.2). The SAS Macro (section A.8.1) runs a data step, so it cannot be called from within another data step, as would usually be desirable. Instead, you would have to generate the data and then merge it into any additional data that was required. On the whole, it might make more sense to treat every application as a special case and hand-code each one, rather than to use the macro, were a built-in function not available.

%macro rtable(outdata=, outvar=, n=1, p=);
data &outdata (keep=&outvar);
/* In the first section of code, we determine the number
   of categories and set up an array to hold the 
   probabilities for each category.  A maximum number of 
   catagories must be provided; here we use 100. */
nprobs = count(&p,",");
lengthprobs = length(&p);
array probs [100] _temporary_;

/* In the next section we parse the parameter string to 
   extract the probabilities for each category. */
lastcomma = 0;
lastprob = 0;
do prob = 1 to nprobs;
 nextcomma = find(&p,",",lastcomma +1);
 probs[prob] = lastprob + 
     substr(&p,lastcomma+1, (nextcomma-lastcomma -1));
 lastcomma = nextcomma;
 lastprob = probs[prob];
 end;
/* Finally, we generate a uniform random variate and
   use it to find the appropriate category. */
do i = 1 to &n;
 __x = uniform(0);
 &outvar = 1;
 do j = 1 to nprobs;
   &outvar = &outvar + (__x gt probs[j]);
       end;
 output;
 end;
run;
%mend rtable;

We've left this code so that an extra comma must be included in the probability parameter string. This could be fixed with an additional line or two in the data step.

We can then run the macro and test the resulting data using proc freq.

%rtable(outdata=test2,outvar= mycat, n = 10000, 
                                   p = ".1,.2,.3,");

proc freq data = test2; tables mycat; run;

                                Cumulative  Cumulative
mycat    Frequency   Percent    Frequency     Percent
------------------------------------------------------
  1        1032       10.32        1032        10.32
  2        2013       20.13        3045        30.45
  3        2956       29.56        6001        60.01
  4        3999       39.99       10000       100.00

R

In contrast, the R code to write a function (section B.5.2) is trivial, using the cut() function as suggested by Douglas Rivers, discussed in section 1.4.10.

rtable1 <- function(cuts,n) {
   output <- cut(runif(n), cuts, labels=FALSE)
   return(output)
}

In fact, the function hardly seems worth writing. To use it, you need to include 0 and 1 as the first and last elements of the cuts vector. Results can be examined using the table() function (section 2.2.1)

> mycat1 <- rtable1(c(0,.1,.3,.6,1),10000)
> summary(factor(mycat1))
   1    2    3    4 
1058 1989 2978 3975 

To add a little bit of programming to it, suppose we want to list the probabilities per category, rather than the cutpoints. We'll omit the final category probability, assuming the input vector sums to less than 1 and has 1 fewer elements than the desired number of categories.

rtable2 <- function(p, n) {
   cuts <- numeric(length(p) + 2)
   cuts[1] <- 0
   cuts[length(p) + 2] <- 1
   for (k in 1:length(p))  {
      cuts[k+1] <- sum(p[1:k])
   }
   output <- cut(runif(n), cuts, labels=FALSE)
   return(output)
}

The results are similar to the simpler function. We can see the results tabulated without explicitly making the result a factor by using the table() function (section 2.2).

> mycat2 <- rtable2(c(.1, .2, .3), 10000)
> table(mycat2)
mycat2
   1    2    3    4 
 961 1986 3058 3995 

Example 7.20: Simulate categorical data

Both SAS and R provide means of simulating categorical data (see section 1.10.4). Alternatively, it is trivial to write code to do this directly. In this entry, we show how to do it once. In a future entry, we'll demonstrate writing a SAS Macro (section A.8.1) and a function in R (section B.5.2) to do it repeatedly.

SAS

data test;
p1 = .1; p2 = .2; p3 = .3;
do i = 1 to 10000;
  x = uniform(0);
  mycat = (x ge 0) + (x gt p1) + (x gt p1 + p2) 
                                + (x gt p1 + p2 + p3);
  output;
  end;
run;

Here the parenthetical logical tests in the mycat = line resolve to 1 if the test is true and 0 otherwise, as discussed in section 1.4.9.
The (x ge 0) makes the categories range from 1 to 4, rather than 0 to 3.

The results can be assessed using proc freq:

proc freq data=test; tables mycat; run;

                                  Cumulative    Cumulative
mycat    Frequency     Percent     Frequency      Percent
----------------------------------------------------------
    1         947        9.47           947         9.47
    2        2061       20.61          3008        30.08
    3        3039       30.39          6047        60.47
    4        3953       39.53         10000       100.00

R

In contrast, the R syntax to get the results is rather dense.

p <- c(.1,.2,.3)
x <- runif(10000)
mycat <- numeric(10000)
for (i in 0:length(p)) { 
    mycat <- mycat + (x >= sum(p[0:i])) 
    }

We can display the results using the summary() function.

summary(factor(mycat))
  1    2    3    4
 990 2047 2978 3985