(This article was first published on

Suppose we wanted to assess the probability P(X=x) for a binomial random variate with n = 10 and with p = .81, .84, ..., .99. This could be helpful, for example, in various game settings. **SAS and R**, and kindly contributed to R-bloggers)In SAS, we ﬁnd the probability that X=x using differences in the CDF calculated via the

`cdf`function (1.10.1). We loop through the various binomial probabilities and values of x using the

`do ... end`structure (section 1.11.1). Finally, we store the probabilities in implicitly named variables via an

`array`statement (section 1.11.5).

**SAS**

data table (drop=j);

array probs [11] prob0 prob1 - prob10;

do p = .81 to .99 by .03;

do j = 0 to 10;

if j eq 1 then probs[j+1] = cdf("BINOMIAL", 0, p, 10);

else probs[j+1] = cdf("BINOMIAL", j, p, 10)

- cdf("BINOMIAL", j-1, p, 10);

end;

output;

end;

run;

The results are printed with two decimal places using the

`format`statement (section 1.2.4). The

`options`statement (section A.4) uses the

`ls`option to narrow the output width to be compatible with Blogger.

options ls=64;

proc print data=table noobs;

var p prob0 prob1 - prob10;

format _numeric_ 3.2;

run;

And the results are:

p

p p p p p p p p p p r

r r r r r r r r r r o

o o o o o o o o o o b

b b b b b b b b b b 1

p 0 1 2 3 4 5 6 7 8 9 0

.81 .00 .00 .00 .00 .00 .02 .08 .19 .30 .29 .12

.84 .00 .00 .00 .00 .00 .01 .05 .15 .29 .33 .17

.87 .00 .00 .00 .00 .00 .00 .03 .10 .25 .37 .25

.90 .00 .00 .00 .00 .00 .00 .01 .06 .19 .39 .35

.93 .00 .00 .00 .00 .00 .00 .00 .02 .12 .36 .48

.96 .00 .00 .00 .00 .00 .00 .00 .01 .05 .28 .66

.99 .00 .00 .00 .00 .00 .00 .00 .00 .00 .09 .90

**R**

In R we start by making a vector of the binomial probabilities, using the

`:`operator (section 1.11.3) to generate a sequence of integers. After creating a matrix (section B.4.4) to hold the table results, we loop (section 1.11.1) through the binomial probabilities, calling the

`dbinom()`function (section 1.1) to ﬁnd the probability that X=x. This calculation is nested within the

`round()`function (section 1.2.5) to reduce the digits displayed. Finally, we glue the vector of binomial probabilities to the results.

p <- .78 + (3 * 1:7)/100

allprobs <- matrix(nrow=length(p), ncol=11)

for (i in 1:length(p)) {

allprobs[i,] <- round(dbinom(0:10, 10, p[i]),2)

}

table <- cbind(p, allprobs)

table

With the result:

p

[1,] 0.81 0 0 0 0 0 0.02 0.08 0.19 0.30 0.29 0.12

[2,] 0.84 0 0 0 0 0 0.01 0.05 0.15 0.29 0.33 0.17

[3,] 0.87 0 0 0 0 0 0.00 0.03 0.10 0.25 0.37 0.25

[4,] 0.90 0 0 0 0 0 0.00 0.01 0.06 0.19 0.39 0.35

[5,] 0.93 0 0 0 0 0 0.00 0.00 0.02 0.12 0.36 0.48

[6,] 0.96 0 0 0 0 0 0.00 0.00 0.01 0.05 0.28 0.66

[7,] 0.99 0 0 0 0 0 0.00 0.00 0.00 0.00 0.09 0.90

As with the example of jittered scatterplots for dichotomous y by continuous x, (Example 7.3, Example 7.4, and Example 7.5) we will revisit this example for improvement in later entries.

To

**leave a comment**for the author, please follow the link and comment on his blog:**SAS and R**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...