Every distribution in R comes with four functions. The name is a letter plus the distribution name:
| Letter | What it does | Answers the question | Example (binomial) |
|---|---|---|---|
d |
Density: probability of exactly one value | P(X = k) | dbinom(6, 20, .16) = P(exactly 6) |
p |
Probability, cumulative: everything up to and including a value | P(X <= k) | pbinom(5, 20, .16) = P(5 or fewer) |
q |
Quantile: the reverse of p; give a
probability, get the value |
“What value has 95% below it?” | qbinom(.95, 20, .16) |
r |
Random: simulate random draws | “Generate fake data” | rbinom(10, 20, .16) = 10 simulated teams |
The distribution names used in this file:
| Distribution | Name in R | Functions | Parameters |
|---|---|---|---|
| Binomial | binom |
dbinom, pbinom, qbinom,
rbinom |
n trials, p probability of success |
| Poisson | pois |
dpois, ppois, qpois,
rpois |
lambda (average count) |
| Hypergeometric | hyper |
dhyper, phyper, qhyper,
rhyper |
m successes, n failures, k drawn |
| Normal | norm |
dnorm, pnorm, qnorm,
rnorm |
mean, sd |
Useful rules:
p functions give the left tail, P(X
<= k), by default. Add lower.tail = FALSE to get the
right tail, P(X > k).pbinom(5, ..., lower.tail = FALSE), using 5, not 6.dnorm gives the height of the curve,
not a probability. Use pnorm for normal probabilities.
Also, <= and < give the same answer.d values gives the same answer as one
p call: sum(dbinom(0:5, 20, .16)) equals
pbinom(5, 20, .16).?dbinom (or ?dpois,
?dhyper, ?dnorm) to open the help page for all
four functions of that distribution.The national 30-day case fatality rate for intracranial hemorrhage (ICH) surgical intervention is 16%. The neurosurgical team performed 20 of these operations last year and experienced 6 deaths. Does sufficient evidence exist to suggest that the neurosurgical practice pattern requires review? In other words, what is the probability that this team would experience something this or more severe, i.e., 6 or more deaths in 20 trials if they were operating at the national rate of .16? Use an evidentiary standard that investigates if the event is rare (routinely defined as less than .05 probability.)
# Brute force: P(X >= 6) = 1 - P(X <= 5).
# Each term is the binomial formula: choose(n, k) * p^k * (1-p)^(n-k),
# where choose(20, k) counts the ways to get k deaths in 20 operations.
1-(choose(20,0)*(.16)^0*(.84)^20 + choose(20,1)*(.16)^1*(.84)^19 + choose(20,2)*(.16)^2*(.84)^18 + choose(20,3)*(.16)^3*(.84)^17 + choose(20,4)*(.16)^4*(.84)^16 + choose(20, 5)*(.16)^5*(.84)^15) #this would be "brute force" method, fully written out## [1] 0.08699685?dbinom # open the help page for the binomial functions (dbinom, pbinom, qbinom, rbinom)
sum(dbinom(6:20, 20, .16)) # dbinom gives P(X = k); add P(X = 6) through P(X = 20) to get P(X >= 6)## [1] 0.086996851-sum(dbinom(0:5, 20, .16)) # same answer via the complement: 1 - [P(X = 0) + ... + P(X = 5)]## [1] 0.08699685dbinom(0:5, 20, .16) # show the six individual probabilities P(X = 0) to P(X = 5)## [1] 0.03059044 0.11653501 0.21087287 0.24099757 0.19509327 0.118913991-pbinom(5, 20, .16) # pbinom gives the cumulative P(X <= 5); 1 minus it gives P(X >= 6)## [1] 0.08699685pbinom(5, 20, .16, lower.tail=FALSE) # lower.tail = FALSE returns P(X > 5) directly, i.e., P(X >= 6)## [1] 0.08699685dbinom(6, 20, .16) # this just illustrates the command, has no bearing on the problem # P(exactly 6 deaths)## [1] 0.05662571The rate of medication errors across the United States is 2 per 1000 orders. (Volume). A sample from your local hospital finds 5 errors in a sample of 2000. What is the probability that this hospital is within the US standard?
?dpois # open the help page for the Poisson functions
# Expected errors in 2000 orders: lambda = 2 per 1000 * 2 = 4.
sum(dpois(5:2000, 4)) # dpois gives P(X = k); add P(X = 5) up to P(X = 2000) to get P(X >= 5)## [1] 0.37116311-ppois(4,4) # complement: 1 - P(X <= 4) = P(X >= 5)## [1] 0.3711631ppois(4, 4, lower.tail=FALSE) # same thing in one step: P(X > 4)## [1] 0.3711631probability = lambda/interval = 2/1000 = .002
# Same question as a binomial: 2000 orders, each with a .002 chance of error.
# With large n and small p, the binomial and Poisson answers nearly match.
pbinom(4, 2000, .002, lower.tail=FALSE) # P(X > 4) = P(X >= 5)## [1] 0.371163Defects occur along the length of cable at an average of six defects per 4000 feet. What is the probability that a 3000-foot cable will have at most two defects?
lambda = 6, t=3000/4000 = .75, lambda_t=6*.75 = 4.5
ppois(2, 4.5) # cumulative P(X <= 2) with lambda = 4.5## [1] 0.1735781sum(dpois(0:2, 4.5)) # same answer: P(X = 0) + P(X = 1) + P(X = 2)## [1] 0.1735781Orthopedic example - 20 patients, 10 will require MRI, what is probability that if doctor sees 5 patients they will send all 5 to MRI
?dhyper # open the help page; arguments are (x, m, n, k)
# x = successes we want, m = successes in the population,
# n = failures in the population, k = number drawn (sampling without replacement)
dhyper(5, 10, 10, 5) # P(all 5 patients seen need an MRI): x = 5, m = 10 MRI, n = 10 no MRI, k = 5 seen## [1] 0.01625387Probability of 4 or more?
sum(dhyper(4:5, 10, 10, 5)) # P(X = 4) + P(X = 5) = P(X >= 4)## [1] 0.1517028phyper(3, 10, 10, 5, lower.tail=FALSE) # same answer: P(X > 3) = P(X >= 4)## [1] 0.151702854 numbers, 6 numbers for jackpot, 3 or more to win something. Total number of successes = 6, number of non-successes = 48, sample 6 numbers
Jackpot = 6 successes out of sample of 6
dhyper(6, 6, 48, 6) # P(match all 6): x = 6, m = 6 winning numbers, n = 48 losing numbers, k = 6 picked## [1] 3.871892e-08Probability of a win (3 or more numbers)
sum(dhyper(3:6, 6, 48, 6)) # P(X = 3) + P(X = 4) + P(X = 5) + P(X = 6)## [1] 0.014059961-phyper(2, 6, 48, 6) # complement: 1 - P(X <= 2)## [1] 0.01405996phyper(2, 6, 48, 6, lower.tail=FALSE) # same thing in one step: P(X > 2)## [1] 0.01405996Mean temperature = 98.6, standard deviation = 1
P(X<96)
?dnorm # open the help page for the normal functions (dnorm, pnorm, qnorm, rnorm)
pnorm(96, 98.6, 1) # pnorm gives the area to the left: P(X < 96) with mean 98.6, sd 1## [1] 0.004661188#via z-score
# Standardize first: z = (x - mean) / sd, then use the standard normal (mean 0, sd 1).
# Here sd = 1, so dividing by sd changes nothing; in general, write (96-98.6)/sd.
pnorm((96-98.6), 0, 1)## [1] 0.004661188P(X>99.5)
1-pnorm(99.5, 98.6, 1) # area to the right: 1 - P(X < 99.5)## [1] 0.1840601pnorm(99.5, 98.6, 1, lower.tail=FALSE) # same thing in one step## [1] 0.1840601P(96<X<98)
pnorm(98, 98.6, 1) - pnorm(96, 98.6, 1) # area between: P(X < 98) - P(X < 96)## [1] 0.2695919P(X>100)
# Area to the right of 100.
# NOTE: the code uses 100.5, which gives P(X > 100.5), not P(X > 100).
# The +0.5 continuity correction applies only when a normal approximates a discrete
# distribution; temperature is continuous. For P(X > 100), use pnorm(100, 98.6, 1, lower.tail=FALSE).
pnorm(100.5, 98.6, 1, lower.tail=FALSE)## [1] 0.02871656