R - Binomial Distribution

R - Binomial Distribution

The binomial distribution is a fundamental discrete probability distribution. It models the number of successes in a fixed number of independent trials, where each trial has only two possible outcomes (e.g., success/failure, heads/tails). For example, we can use it to calculate the probability of getting exactly 7 heads in 10 coin tosses.

R provides a suite of four functions for working with the binomial distribution. These functions are powerful tools for both calculation and simulation.

# Binomial distribution functions in R
dbinom(x, size, prob)   # Density function (PDF) - P(X = x)
pbinom(q, size, prob)   # Cumulative distribution function (CDF) - P(X ≤ q)
qbinom(p, size, prob)   # Quantile function (inverse CDF)
rbinom(n, size, prob)   # Random number generation

Here’s a breakdown of the parameters:

• x, q: Vector of quantiles (the number of successes).
• p: Vector of probabilities.
• n: Number of random values to generate.
• size: The number of trials (e.g., 10 coin tosses).
• prob: The probability of success on a single trial (e.g., 0.5 for a fair coin).

dbinom() - Probability Mass Function

This function calculates the probability of observing exactly x successes in size trials. Let’s visualize the probability distribution for a scenario with 50 trials and a success probability of 0.5 (like 50 coin flips).

For modern, publication-quality graphics, we’ll use the ggplot2 package from the tidyverse.

# Ensure you have tidyverse installed: install.packages("tidyverse")
library(ggplot2)
library(dplyr)

# Define parameters
num_trials <- 50
prob_success <- 0.5

# Create a data frame for plotting
binom_data <- tibble(
  successes = 0:num_trials,
  probability = dbinom(x = successes, size = num_trials, prob = prob_success)
)

# Create the plot using ggplot2
binom_plot <- ggplot(binom_data, aes(x = successes, y = probability)) +
  geom_col(fill = "#0072B2", alpha = 0.8) + # Use columns for discrete distribution
  labs(
    title = "Binomial Probability Distribution",
    subtitle = paste("n = ", num_trials, ", p = ", prob_success),
    x = "Number of Successes",
    y = "Probability"
  ) +
  theme_minimal() # Use a clean theme

# Display the plot
print(binom_plot)

# To save the plot to a file:
# ggsave("binomial_distribution.png", plot = binom_plot, width = 8, height = 5)

When you execute this code, it produces a bar chart. The x-axis shows the number of successes (from 0 to 50), and the y-axis shows the probability of each outcome. The chart will be bell-shaped, peaking at 25 successes, which is the most likely outcome.

pbinom() - Cumulative Probability

This function calculates the cumulative probability of an event, i.e., the probability of getting q or fewer successes. This is useful for answering questions like “What is the probability of getting 26 or fewer heads in 51 coin tosses?”

# Calculate the probability of getting 26 or fewer heads from 51 tosses of a fair coin.
prob_leq_26 <- pbinom(q = 26, size = 51, prob = 0.5)

print(prob_leq_26)

When we execute the above code, it produces the following result:

[1] 0.610116

qbinom() - The Quantile Function

This is the inverse of pbinom(). It takes a probability p and returns the number of successes x such that pbinom(x, ...) is equal to p. It answers questions like: “How many heads would mark the 25th percentile of outcomes when a coin is tossed 51 times?”

# Find the number of heads below which 25% of the outcomes fall
# when a coin is tossed 51 times.
quantile_val <- qbinom(p = 0.25, size = 51, prob = 0.5)

print(quantile_val)

When we execute the above code, it produces the following result:

[1] 23

rbinom() - Random Number Generation

This function is used for simulation. It generates n random values from a binomial distribution with specified parameters. This is incredibly useful for simulating experiments.

# Generate 8 random outcomes from an experiment with 150 trials
# and a success probability of 0.4.
set.seed(42) # Set a seed for reproducibility
random_outcomes <- rbinom(n = 8, size = 150, prob = 0.4)

print(random_outcomes)

Because the output is random, your results will differ unless you use set.seed(). With the seed set to 42, the result is:

[1] 61 67 56 61 64 61 58 59