Random Variables: Understand discrete and continuous random variables and their probability distributions
Probability Distributions: Master key continuous distributions (normal, exponential, Weibull, gamma, beta, lognormal) and discrete distributions (binomial, Poisson)
Distribution Properties: Calculate means, variances, and probabilities using probability density/mass functions and cumulative distribution functions
Normal Distribution: Apply standardization techniques and use normal approximations to other distributions
Engineering Applications: Model real-world engineering phenomena using appropriate probability distributions
Multiple Variables: Analyze joint distributions, independence, and functions of random variables including error propagation
After careful study of this chapter, you should be able to do the following:
Determine probabilities for discrete random variables from probability mass functions and for continuous random variables from probability density functions, and use cumulative distribution functions in both cases.
Calculate means and variances for discrete and continuous random variables.
Understand the assumptions for each of the probability distributions presented.
Select an appropriate probability distribution to calculate probabilities in specific applications.
Use the table (or software) for the cumulative distribution function of a standard normal distribution to calculate probabilities.
Approximate probabilities for binomial and Poisson distributions.
Interpret and calculate covariances and correlations between random variables.
Calculate means and variances for linear combinations of random variables.
Approximate means and variances for general functions of several random variables using error propagation.
Understand statistics and the central limit theorem.
A random experiment is an experiment that can result in different outcomes, even though it is repeated in the same manner every time. The collection of all possible outcomes of a random experiment is called the sample space of the experiment.
Key Concepts:
Figure 1 illustrates the relationship between a mathematical model and the physical system it represents. While models like Newton’s laws are not perfect abstractions, they are useful for analyzing and approximating system performance. Once validated with measurements, models can help understand, describe, and predict a system’s response to inputs.
Figure 2 illustrates a model where uncontrollable variables (noise) interact with controllable variables to produce system outputs. Due to the noise, identical controllable settings yield varying outputs upon measurement.
For measuring current in a copper wire, Ohm’s law may suffice as a model. However, if variations are significant, the model may need to account for them (see Figure 3).
In designing a voice communication system, a model is required for call frequency and duration. Even if calls occur and last precisely 5 minutes on average, variations in timing or duration can lead to call blocking, necessitating multiple lines (see Figure 4).
A random variable is a function that assigns a real number to each outcome in the sample space of a random experiment. It provides a numerical description of the outcome of a random experiment.
Definition: A random variable is a numerical variable whose measured value can change from one replicate of the experiment to another.
Types of Random Variables:
The distinction between discrete and continuous random variables is important because different methods are used to work with each type.
Probability is used to quantify likelihood or chance and represents risk or uncertainty in engineering applications. It can be interpreted as our degree of belief or relative frequency.
Key Properties:
For example, if X denotes part length, the probability statement can be written as: P(X \in [10.8, 11.2]) = 0.25 \quad \text{or} \quad P(10.8 \leq X \leq 11.2) = 0.25
Both equations state that the probability that the random variable X assumes a value in [10.8, 11.2] is 0.25.
Given a set E, the complement of E is the set of elements that are not in E. The complement is denoted as E' or E^c.
Property: P(E^c) = 1 - P(E)
The sets E_1, E_2, \ldots, E_k are mutually exclusive if the intersection of any pair is empty. That is, each element is in one and only one of the sets E_1, E_2, \ldots, E_k.
Property: If events are mutually exclusive, their probabilities add.
\small P(X \in E_1 \cup E_2 \cup \ldots \cup E_k) = P(X \in E_1) + P(X \in E_2) + \cdots + P(X \in E_k)
An event is a subset of the sample space. Sometimes, experimental results are classified into categories rather than measured numerically.
Examples:
Events provide a framework for probability calculations in both discrete and continuous settings.
The probability distribution or simply distribution of a random variable X is a description of the set of probabilities associated with the possible values for X.
The probability density function (PDF) f(x) of a continuous random variable X is used to determine probabilities as follows:
P(a < X < b) = \int_{a}^{b} f(x)\,dx
Properties of the PDF:
Important Note: If X is a continuous random variable, for any x_1 and x_2:
\small P(x_1 \leq X \leq x_2) = P(x_1 < X \leq x_2) = P(x_1 \leq X < x_2) = P(x_1 < X < x_2)
This is because P(X = c) = 0 for any specific value c when X is continuous.
Let the continuous random variable X denote the distance in micrometers from the start of a track on a magnetic disk until the first flaw. Historical data show that the distribution of X can be modeled by the probability density function:
f(x) = \frac{1}{2000} e^{-x/2000}, \quad x \geq 0
Problem 1: For what proportion of disks is the distance to the first flaw greater than 1000 micrometers?
Solution: \begin{aligned} P(X > 1000) &= \int_{1000}^{\infty} f(x)\,dx \\ &= \int_{1000}^{\infty} \frac{1}{2000} e^{-x/2000} \,dx \\ &= \left[ -e^{-x/2000} \right]_{1000}^{\infty} \\ &= 0 - (-e^{-1000/2000}) \\ &= e^{-1/2} \\ &\approx 0.6065 \end{aligned}
Problem 2: What proportion of disks has the flaw distance between 1000 and 2000 micrometers?
Solution: \begin{aligned} P(1000 \leq X \leq 2000) &= \int_{1000}^{2000} f(x)\,dx \\ &= \left[ -e^{-x/2000} \right]_{1000}^{2000} \\ &= e^{-1/2} - e^{-1} \\ &= 0.6065 - 0.3679 \\ &= 0.2386 \end{aligned}
[1] "Disk Flaw Distance Analysis:" lambda mean_param prob_greater_1000 prob_between_1000_2000
<num> <num> <num> <num>
1: 5e-04 2000 0.6065307 0.2386512
prob_greater_1000_manual prob_between_manual
<num> <num>
1: 0.6065307 0.2386512# Example 3-1: Exponential Distribution for Disk Flaw Distance
library(fastverse)
library(tidyverse)
library(fitdistrplus)
library(scales)
# Define the exponential PDF for disk flaw example
disk_data <- data.table(
lambda = 1 / 2000,
mean_param = 2000
) %>%
fmutate(
# Calculate probabilities
prob_greater_1000 = pexp(q = 1000, rate = lambda, lower.tail = FALSE),
prob_between_1000_2000 = pexp(q = 2000, rate = lambda, lower.tail = TRUE) -
pexp(q = 1000, rate = lambda, lower.tail = TRUE),
# Verify using manual integration
prob_greater_1000_manual = exp(-1000 / mean_param),
prob_between_manual = exp(-1000 / mean_param) - exp(-2000 / mean_param)
)
print("Disk Flaw Distance Analysis:")
print(disk_data)
# Create visualization
x_vals <- seq(0, 8000, by = 50)
pdf_vals <- dexp(x_vals, rate = 1 / 2000)
fwrite(data.table(x = x_vals, pdf = pdf_vals), "data/exponential_pdf.csv")The cumulative distribution function (CDF) of a continuous random variable X with probability density function f(x) is:
F(x) = P(X \leq x) = \int_{-\infty}^{x} f(u)\,du
for -\infty < x < \infty.
Key Relationships:
Continuing with the disk flaw example, the CDF is:
\begin{aligned} F(x) &= P(X \leq x) = \int_{0}^{x} \frac{1}{2000} e^{-u/2000}\,du \\ &= \left[ -e^{-u/2000} \right]_{0}^{x} \\ &= 1 - e^{-x/2000} \end{aligned}
for x \geq 0, and F(x) = 0 for x < 0.
Verification: We can check that P(X > 1000) = 1 - F(1000) = 1 - (1 - e^{-1/2}) = e^{-1/2} = 0.6065
[1] "CDF Verification:" x cdf_value cdf_manual prob_greater_1000_cdf prob_between_1000_2000_cdf
<num> <num> <num> <num> <num>
1: 1000 0.3934693 0.3934693 0.6065307 0.2386512
2: 2000 0.6321206 0.6321206 0.6065307 0.2386512# Example 3-2: CDF for Disk Flaw Distance
cdf_data <- data.table(
x = x_vals
) %>%
fmutate(
cdf_value = pexp(x, rate = 1 / 2000),
cdf_manual = 1 - exp(-x / 2000),
# Verify probabilities using CDF
prob_greater_1000_cdf = 1 - pexp(1000, rate = 1 / 2000),
prob_between_1000_2000_cdf = pexp(2000, rate = 1 / 2000) - pexp(1000, rate = 1 / 2000)
)
print("CDF Verification:")
print(cdf_data[x %in% c(1000, 2000)])Suppose X is a continuous random variable with PDF f(x).
Mean (Expected Value): The mean or expected value of X, denoted as \mu or E(X), is:
\mu = E(X) = \int_{-\infty}^{\infty} x f(x)\,dx
Variance: The variance of X, denoted as V(X) or \sigma^2, is:
\begin{aligned} \sigma^2 = V(X) &= \int_{-\infty}^{\infty} (x-\mu)^2 f(x)\,dx \\ &= E(X^2) - [E(X)]^2 \end{aligned}
Standard Deviation: The standard deviation of X is \sigma = \sqrt{V(X)}.
Interpretation:
For the disk flaw distance example with f(x) = \frac{1}{2000} e^{-x/2000} for x \geq 0:
Mean Calculation: \begin{aligned} E(X) &= \int_{0}^{\infty} x \cdot \frac{1}{2000} e^{-x/2000}\,dx \end{aligned}
Using integration by parts with u = x and dv = \frac{1}{2000} e^{-x/2000} dx:
\begin{aligned} E(X) &= \left[ -x e^{-x/2000} \right]_{0}^{\infty} + \int_{0}^{\infty} e^{-x/2000}\,dx \\ &= 0 + 2000 \left[ -e^{-x/2000} \right]_{0}^{\infty} \\ &= 2000(0 - (-1)) = 2000 \end{aligned}
Variance Calculation: For this exponential distribution, V(X) = (2000)^2 = 4,000,000.
Interpretation: On average, the first flaw occurs at 2000 micrometers, with substantial variability (standard deviation = 2000 micrometers).
[1] "Exponential Distribution Statistics:" rate scale theoretical_mean theoretical_variance theoretical_sd mean_manual
<num> <num> <num> <num> <num> <num>
1: 5e-04 2000 2000 4e+06 2000 2000
variance_manual
<num>
1: 4e+06# Example 3-3: Mean and Variance Calculation
exponential_stats <- data.table(
rate = 1 / 2000,
scale = 2000
) %>%
fmutate(
theoretical_mean = 1 / rate,
theoretical_variance = 1 / rate^2,
theoretical_sd = sqrt(theoretical_variance),
# Verify mean calculation manually
mean_manual = scale,
variance_manual = scale^2
)
print("Exponential Distribution Statistics:")
print(exponential_stats)The normal distribution is the most widely used model for continuous random variables. It’s also called the Gaussian distribution.
Probability Density Function: A random variable X has a normal distribution with parameters \mu and \sigma if:
f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right\}
for -\infty < x < \infty, where -\infty < \mu < \infty and \sigma > 0.
Parameters:
Notation: X \sim N(\mu, \sigma^2)
Key Properties:
Assume that current measurements in a wire follow a normal distribution with mean \mu = 10 milliamperes and variance \sigma^2 = 4 (so \sigma = 2) milliamperes^2. What is the probability that a measurement exceeds 13 milliamperes?
Solution: We need P(X > 13) where X \sim N(10, 4).
This probability corresponds to the area under the normal curve to the right of 13, as shown in Figure 8.
Since there’s no closed-form expression for the normal integral, we use standardization or statistical software.
For any normal distribution:
\begin{aligned} P(\mu - \sigma < X < \mu + \sigma) &= 0.6827 \\ P(\mu - 2\sigma < X < \mu + 2\sigma) &= 0.9545 \\ P(\mu - 3\sigma < X < \mu + 3\sigma) &= 0.9973 \end{aligned}
[1] "Normal Distribution Empirical Rule:" mu sigma prob_1sigma prob_2sigma prob_3sigma theoretical_1sigma
<num> <num> <num> <num> <num> <num>
1: 0 1 0.6826895 0.9544997 0.9973002 0.6827
theoretical_2sigma theoretical_3sigma
<num> <num>
1: 0.9545 0.9973# Normal Distribution Properties and Empirical Rule
normal_props <- data.table(
mu = 0,
sigma = 1
) %>%
fmutate(
# Empirical rule calculations
prob_1sigma = pnorm(mu + sigma) - pnorm(mu - sigma),
prob_2sigma = pnorm(mu + 2 * sigma) - pnorm(mu - 2 * sigma),
prob_3sigma = pnorm(mu + 3 * sigma) - pnorm(mu - 3 * sigma),
# Theoretical values
theoretical_1sigma = 0.6827,
theoretical_2sigma = 0.9545,
theoretical_3sigma = 0.9973
)
print("Normal Distribution Empirical Rule:")
print(normal_props)A normal random variable with \mu = 0 and \sigma^2 = 1 is called a standard normal random variable, denoted as Z.
Standard Normal CDF: \Phi(z) = P(Z \leq z)
Standardizing: If X \sim N(\mu, \sigma^2), then:
Z = \frac{X - \mu}{\sigma} \sim N(0, 1)
Key Property: P(X \leq x) = P\left(Z \leq \frac{x - \mu}{\sigma}\right) = \Phi\left(\frac{x - \mu}{\sigma}\right)
Problem 1: P(Z > 1.26)
Solution: P(Z > 1.26) = 1 - P(Z \leq 1.26) = 1 - 0.89617 = 0.10383
Problem 2: P(Z < -0.86)
Solution: P(Z < -0.86) = 0.19489
Problem 3: P(-1.25 < Z < 0.37)
Solution: P(-1.25 < Z < 0.37) = P(Z < 0.37) - P(Z < -1.25) = 0.64431 - 0.10565 = 0.53866
Problem 4: Find z such that P(Z > z) = 0.05
Solution: We need P(Z \leq z) = 0.95. From the table, z = 1.645.
[1] "Standard Normal Calculations:" calculation result
<char> <num>
1: P(Z > 1.26) 1.038347e-01
2: P(Z < -0.86) 1.948945e-01
3: P(Z > -1.37) 9.146565e-01
4: P(-1.25 < Z < 0.37) 5.386590e-01
5: P(Z <= -4.6) 2.112455e-06
6: z for P(Z > z) = 0.05 1.644854e+00
7: z for P(-z < Z < z) = 0.99 2.575829e+00# Example 3-5: Standard Normal Calculations
standard_normal_calcs <- data.table(
calculation = c(
"P(Z > 1.26)", "P(Z < -0.86)", "P(Z > -1.37)",
"P(-1.25 < Z < 0.37)", "P(Z <= -4.6)", "z for P(Z > z) = 0.05",
"z for P(-z < Z < z) = 0.99"
)
) %>%
fmutate(
result = case_when(
calculation == "P(Z > 1.26)" ~ pnorm(1.26, lower.tail = FALSE),
calculation == "P(Z < -0.86)" ~ pnorm(-0.86),
calculation == "P(Z > -1.37)" ~ pnorm(-1.37, lower.tail = FALSE),
calculation == "P(-1.25 < Z < 0.37)" ~ pnorm(0.37) - pnorm(-1.25),
calculation == "P(Z <= -4.6)" ~ pnorm(-4.6),
calculation == "z for P(Z > z) = 0.05" ~ qnorm(0.05, lower.tail = FALSE),
calculation == "z for P(-z < Z < z) = 0.99" ~ qnorm(0.995),
TRUE ~ NA_real_
)
)
print("Standard Normal Calculations:")
print(standard_normal_calcs)The diameter of a shaft is normally distributed with mean 0.2508 inch and standard deviation 0.0005 inch. The specifications are 0.2500 \pm 0.0015 inch. What proportion of shafts conforms to specifications?
Solution: We need P(0.2485 < X < 0.2515) where X \sim N(0.2508, 0.0005^2).
Standardizing: \small \begin{aligned} P(0.2485 < X < 0.2515) &= P\left(\frac{0.2485 - 0.2508}{0.0005} < Z < \frac{0.2515 - 0.2508}{0.0005}\right) \\ &= P(-4.6 < Z < 1.4) \\ &= P(Z < 1.4) - P(Z < -4.6) \\ &= 0.91924 - 0.00000 = 0.91924 \end{aligned}
Interpretation: About 91.9% of shafts meet specifications.
[1] "Shaft Diameter Quality Control:" mu sigma lower_spec upper_spec z_lower z_upper prob_meets_specs
<num> <num> <num> <num> <num> <num> <num>
1: 0.2508 5e-04 0.2485 0.2515 -4.6 1.4 0.9192412
prob_direct
<num>
1: 0.9192412# Example 3-6: Shaft Diameter Quality Control
shaft_data <- data.table(
mu = 0.2508,
sigma = 0.0005,
lower_spec = 0.2485,
upper_spec = 0.2515
) %>%
fmutate(
# Standardize the specification limits
z_lower = (lower_spec - mu) / sigma,
z_upper = (upper_spec - mu) / sigma,
# Calculate probability
prob_meets_specs = pnorm(z_upper) - pnorm(z_lower),
# Direct calculation
prob_direct = pnorm(upper_spec, mean = mu, sd = sigma) -
pnorm(lower_spec, mean = mu, sd = sigma)
)
print("Shaft Diameter Quality Control:")
print(shaft_data)If W \sim N(\theta, \omega^2), then X = \exp(W) has a lognormal distribution.
Probability Density Function: f(x) = \frac{1}{x\sqrt{2\pi\omega^2}} \exp\left\{-\frac{1}{2}\left(\frac{\ln(x) - \theta}{\omega}\right)^2\right\}
for 0 < x < \infty.
Parameters:
Mean and Variance:
\begin{aligned} E(X) &= \exp\left(\theta + \frac{\omega^2}{2}\right) \\ V(X) &= \exp(2\theta + \omega^2)\left[\exp(\omega^2) - 1\right] \end{aligned}
Applications: Lifetimes, financial data, environmental measurements
The lifetime of a semiconductor laser has a lognormal distribution with \theta = 10 and \omega = 1.5 hours.
Problem 1: What is the probability the lifetime exceeds 10,000 hours?
Solution: Let W \sim N(10, 1.5^2), so X = \exp(W).
\begin{aligned} P(X > 10,000) &= P(\exp(W) > 10,000) \\ &= P(W > \ln(10,000)) \\ &= P\left(\frac{W - 10}{1.5} > \frac{\ln(10,000) - 10}{1.5}\right) \\ &= P\left(Z > \frac{9.2103 - 10}{1.5}\right) \\ &= P(Z > -0.5264) = 0.7009 \end{aligned}
Problem 2: What lifetime is exceeded by 99% of lasers?
Solution: We need x such that P(X > x) = 0.99, or P(X \leq x) = 0.01.
From P(\ln(X) \leq \ln(x)) = 0.01 and \ln(X) \sim N(10, 1.5^2):
P\left(\frac{\ln(X) - 10}{1.5} \leq \frac{\ln(x) - 10}{1.5}\right) = 0.01
Since P(Z \leq -2.326) = 0.01: \frac{\ln(x) - 10}{1.5} = -2.326 x = \exp(10 - 1.5 \times 2.326) = \exp(6.511) = 673.6 \text{ hours}
[1] "Lognormal Distribution Analysis:" theta omega test_value prob_exceeds_10000 percentile_99 theoretical_mean
<num> <num> <num> <num> <num> <num>
1: 10 1.5 10000 0.7007086 672.1478 67846.29
theoretical_variance theoretical_sd
<num> <num>
1: 39070059887 197661.5# Example 3-7: Semiconductor Laser Lifetime (Lognormal)
lognormal_data <- data.table(
theta = 10,
omega = 1.5,
test_value = 10000
) %>%
fmutate(
# Probability calculations
prob_exceeds_10000 = plnorm(test_value, meanlog = theta, sdlog = omega, lower.tail = FALSE),
# Find 99th percentile (exceeded by 99%)
percentile_99 = qlnorm(0.01, meanlog = theta, sdlog = omega),
# Mean and variance
theoretical_mean = exp(theta + omega^2 / 2),
theoretical_variance = exp(2 * theta + omega^2) * (exp(omega^2) - 1),
theoretical_sd = sqrt(theoretical_variance)
)
print("Lognormal Distribution Analysis:")
print(lognormal_data)The random variable X has a gamma distribution with shape parameter r > 0 and rate parameter \lambda > 0 if:
Probability Density Function:
f(x) = \frac{\lambda^r x^{r-1} \exp(-\lambda x)}{\Gamma(r)}, \quad x > 0
where \Gamma(r) = \int_0^{\infty} t^{r-1} e^{-t} dt is the gamma function.
Properties of Gamma Function:
Mean and Variance: \begin{aligned} E(X) &= \frac{r}{\lambda} \\ V(X) &= \frac{r}{\lambda^2} \end{aligned}
Special Cases:
The Weibull distribution is widely used in reliability engineering and survival analysis.
Probability Density Function: f(x) = \frac{\beta}{\delta} \left(\frac{x}{\delta}\right)^{\beta-1} \exp\left\{-\left(\frac{x}{\delta}\right)^{\beta}\right\}, \quad x > 0
Cumulative Distribution Function: F(x) = 1 - \exp\left\{-\left(\frac{x}{\delta}\right)^{\beta}\right\}
Parameters:
Mean and Variance: \begin{aligned} E(X) &= \delta \Gamma\left(1 + \frac{1}{\beta}\right) \\ V(X) &= \delta^2 \left[\Gamma\left(1 + \frac{2}{\beta}\right) - \left\{\Gamma\left(1 + \frac{1}{\beta}\right)\right\}^2\right] \end{aligned}
Special Cases:
The time to failure (in hours) of a bearing follows a Weibull distribution with \beta = 0.5 and \delta = 5000 hours.
Problem 1: Determine the mean time until failure.
Solution: \begin{aligned} E(X) &= \delta \Gamma\left(1 + \frac{1}{\beta}\right) \\ &= 5000 \times \Gamma(1 + 2) \\ &= 5000 \times \Gamma(3) \\ &= 5000 \times 2! = 10,000 \text{ hours} \end{aligned}
Problem 2: Determine the probability that a bearing lasts at least 6000 hours.
Solution: \begin{aligned} P(X > 6000) &= 1 - F(6000) \\ &= \exp\left\{-\left(\frac{6000}{5000}\right)^{0.5}\right\} \\ &= \exp(-1.2^{0.5}) \\ &= \exp(-1.095) = 0.334 \end{aligned}
Interpretation: Only 33.4% of bearings last at least 6000 hours.
[1] "Weibull Distribution Analysis:" shape scale test_time mean_failure_time prob_lasts_6000 variance sd
<num> <num> <num> <num> <num> <num> <num>
1: 0.5 5000 6000 10000 0.3343907 5e+08 22360.68# Example 3-8: Bearing Failure Time (Weibull)
weibull_data <- data.table(
shape = 0.5,
scale = 5000,
test_time = 6000
) %>%
fmutate(
# Mean time until failure
mean_failure_time = scale * gamma(1 + 1 / shape),
# Probability of lasting at least 6000 hours
prob_lasts_6000 = pweibull(test_time, shape = shape, scale = scale, lower.tail = FALSE),
# Variance calculation
variance = scale^2 * (gamma(1 + 2 / shape) - (gamma(1 + 1 / shape))^2),
sd = sqrt(variance)
)
print("Weibull Distribution Analysis:")
print(weibull_data)The beta distribution is defined on the interval [0, 1] and is useful for modeling proportions and percentages.
Probability Density Function:
f(x) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1} (1-x)^{\beta-1}, \quad 0 \leq x \leq 1
Parameters:
Mean and Variance:
\begin{aligned} E(X) &= \frac{\alpha}{\alpha + \beta} \\ V(X) &= \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \end{aligned}
Special Cases:
Applications: Quality control, project completion times, proportions
The proportion of maximum allowed time to complete a task follows a beta distribution with \alpha = 2.5 and \beta = 1.
Problem 1: What is the probability that the proportion exceeds 0.7?
Solution: \begin{aligned} P(X > 0.7) &= \int_{0.7}^{1} \frac{\Gamma(2.5 + 1)}{\Gamma(2.5)\Gamma(1)} x^{2.5-1} (1-x)^{1-1} dx \\ &= \int_{0.7}^{1} \frac{\Gamma(3.5)}{\Gamma(2.5)} x^{1.5} dx \\ &= \frac{2.5}{1} \int_{0.7}^{1} x^{1.5} dx \\ &= 2.5 \left[\frac{x^{2.5}}{2.5}\right]_{0.7}^{1} \\ &= (1^{2.5} - 0.7^{2.5}) = 1 - 0.7^{2.5} = 0.59 \end{aligned}
Problem 2: Calculate the mean and variance.
Solution: \begin{aligned} E(X) &= \frac{\alpha}{\alpha + \beta} = \frac{2.5}{2.5 + 1} = \frac{2.5}{3.5} = 0.714 \\ V(X) &= \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)} \\ &= \frac{2.5 \times 1}{(3.5)^2 \times 4.5} = \frac{2.5}{55.125} = 0.045 \end{aligned}
[1] "Beta Distribution Analysis:" alpha beta_param test_value prob_exceeds_07 mean_beta variance_beta
<num> <num> <num> <num> <num> <num>
1: 2.5 1 0.7 0.5900366 0.7142857 0.04535147
sd_beta
<num>
1: 0.2129589# Example 3-9: Project Completion Time (Beta)
beta_data <- data.table(
alpha = 2.5,
beta_param = 1,
test_value = 0.7
) %>%
fmutate(
# Probability exceeds 0.7
prob_exceeds_07 = pbeta(test_value, shape1 = alpha, shape2 = beta_param, lower.tail = FALSE),
# Mean and variance
mean_beta = alpha / (alpha + beta_param),
variance_beta = (alpha * beta_param) / ((alpha + beta_param)^2 * (alpha + beta_param + 1)),
sd_beta = sqrt(variance_beta)
)
print("Beta Distribution Analysis:")
print(beta_data)A normal probability plot is a graphical method for determining whether sample data conform to a hypothesized normal distribution. The procedure is:
Order the sample data from smallest to largest: x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}
Plot the pairs (x_{(i)}, \Phi^{-1}((i-0.5)/n)) for i = 1, 2, \ldots, n
If the data are from a normal distribution, the plotted points will approximately follow a straight line
Interpretation:
Advantages:
Consider battery life data (in hours): 176, 191, 214, 220, 205, 192, 201, 190, 183, 185.
The normal probability plot helps assess whether these data could reasonably come from a normal distribution.
[1] "Battery Life Data for Normal Probability Plot:" life sample_mean sample_sd sample_size rank theoretical_quantile
<num> <num> <num> <int> <num> <num>
1: 176 195.7 14.03211 10 1 -1.6448536
2: 183 195.7 14.03211 10 2 -1.0364334
3: 185 195.7 14.03211 10 3 -0.6744898
4: 190 195.7 14.03211 10 4 -0.3853205
5: 191 195.7 14.03211 10 5 -0.1256613
6: 192 195.7 14.03211 10 6 0.1256613
7: 201 195.7 14.03211 10 7 0.3853205
8: 205 195.7 14.03211 10 8 0.6744898
9: 214 195.7 14.03211 10 9 1.0364334
10: 220 195.7 14.03211 10 10 1.6448536
fitted_value
<num>
1: 172.6192
2: 181.1567
3: 186.2355
4: 190.2931
5: 193.9367
6: 197.4633
7: 201.1069
8: 205.1645
9: 210.2433
10: 218.7808[1] "Shapiro-Wilk p-value: 0.7173"# Example 3-10: Normal Probability Plot
battery_life <- data.table(
life = c(176, 191, 214, 220, 205, 192, 201, 190, 183, 185)
) %>%
fmutate(
sample_mean = fmean(life),
sample_sd = fsd(life),
sample_size = fnobs(life)
)
# Create normal probability plot data
battery_sorted <- battery_life %>%
fmutate(rank = frank(life, ties.method = "average")) %>%
fmutate(
theoretical_quantile = qnorm((rank - 0.5) / sample_size),
fitted_value = sample_mean + sample_sd * theoretical_quantile
) %>%
roworder(life)
print("Battery Life Data for Normal Probability Plot:")
print(battery_sorted)
# Test for normality using Shapiro-Wilk
shapiro_test <- shapiro.test(battery_life$life)
print(paste("Shapiro-Wilk p-value:", round(shapiro_test$p.value, 4)))
fwrite(battery_sorted, "data/battery_normal_plot.csv")Interpretation: If the points roughly follow a straight line, the normal distribution is a reasonable model for the data.
The same probability plotting technique can be used with any distribution:
Lognormal Probability Plot:
Weibull Probability Plot:
Exponential Probability Plot:
General Approach:
Consider failure time data that might follow a lognormal distribution. The lognormal probability plot helps determine if this distribution is appropriate.
[1] "Lognormal Distribution Fit:" time ln_time sample_size meanlog_est sdlog_est
<num> <num> <int> <num> <num>
1: 81 4.394449 50 4.966551 0.4756993[1] "KS test p-value for lognormal: 0.8849"# Example 3-11: Lognormal Probability Plot
# Simulate some failure time data that might be lognormal
set.seed(42)
failure_times <- data.table(
time = c(
81, 249, 117, 227, 134, 98, 135, 149, 225, 59,
291, 223, 127, 185, 181, 101, 205, 115, 240, 151,
98, 80, 198, 161, 240, 118, 177, 342, 197, 146,
158, 82, 83, 98, 104, 197, 64, 34, 65, 100,
139, 137, 342, 144, 215, 249, 149, 185, 151, 200
)
) %>%
fmutate(
ln_time = log(time),
sample_size = fnobs(time)
)
# Fit lognormal distribution
lognormal_fit <- failure_times %>%
fmutate(
meanlog_est = fmean(ln_time),
sdlog_est = fsd(ln_time)
)
print("Lognormal Distribution Fit:")
print(lognormal_fit[1])
# Test goodness of fit
ks_test_lognormal <- ks.test(failure_times$time, "plnorm",
meanlog = lognormal_fit$meanlog_est[1],
sdlog = lognormal_fit$sdlog_est[1]
)
print(paste("KS test p-value for lognormal:", round(ks_test_lognormal$p.value, 4)))
fwrite(failure_times, "data/failure_times.csv")A discrete random variable can take on only a finite or countably infinite number of distinct values. Common examples include:
Key Characteristics:
For a discrete random variable X with possible values x_1, x_2, \ldots, x_n, the probability mass function is:
f(x_i) = P(X = x_i)
Properties:
Notation: Sometimes written as p(x) or P(X = x)
Let X equal the number of bits in error in the next 4 bits transmitted. Suppose the probabilities are:
This completely specifies the probability distribution of X.
[1] "Transmission Error Distribution:" errors probability total_prob cumulative_prob
<int> <num> <num> <num>
1: 0 0.6561 1 0.6561
2: 1 0.2916 1 0.9477
3: 2 0.0486 1 0.9963
4: 3 0.0036 1 0.9999
5: 4 0.0001 1 1.0000# Example 3-12: Digital Transmission Errors (PMF)
transmission_data <- data.table(
errors = 0:4,
probability = c(0.6561, 0.2916, 0.0486, 0.0036, 0.0001)
) %>%
fmutate(
# Verify probabilities sum to 1
total_prob = fsum(probability),
# Calculate cumulative probabilities
cumulative_prob = cumsum(probability)
)
print("Transmission Error Distribution:")
print(transmission_data)
fwrite(transmission_data, "data/transmission_errors.csv")The cumulative distribution function of a discrete random variable X is:
F(x) = P(X \leq x) = \sum_{x_i \leq x} f(x_i)
Properties:
Key Relationship: P(a < X \leq b) = F(b) - F(a)
Continuing the previous example:
Note: F(1.5) = P(X \leq 1.5) = P(X \leq 1) = 0.9477 since X cannot equal 1.5.
[1] "CDF for Transmission Errors:" errors probability cumulative_prob prob_statement
<int> <num> <num> <char>
1: 0 0.6561 0.6561 P(X <= 0) = 0.6561
2: 1 0.2916 0.9477 P(X <= 1) = 0.9477
3: 2 0.0486 0.9963 P(X <= 2) = 0.9963
4: 3 0.0036 0.9999 P(X <= 3) = 0.9999
5: 4 0.0001 1.0000 P(X <= 4) = 1[1] "CDF at Non-Integer Values:" x cdf_value interpretation
<num> <num> <char>
1: 1.5 0.9477 F(1.5) = 0.9477
2: 2.3 0.9963 F(2.3) = 0.9963
3: 3.7 0.9999 F(3.7) = 0.9999# Example 3-13: CDF for Transmission Errors
cdf_transmission <- transmission_data %>%
fselect(errors, probability, cumulative_prob) %>%
fmutate(
prob_statement = paste0("P(X <= ", errors, ") = ", round(cumulative_prob, 4))
)
print("CDF for Transmission Errors:")
print(cdf_transmission)
# Example of intermediate values
intermediate_vals <- data.table(
x = c(1.5, 2.3, 3.7),
cdf_value = c(0.9477, 0.9963, 0.9999)
) %>%
fmutate(
interpretation = paste0("F(", x, ") = ", cdf_value)
)
print("CDF at Non-Integer Values:")
print(intermediate_vals)Let X have possible values x_1, x_2, \ldots, x_n with PMF f(x).
Mean (Expected Value):
\mu = E(X) = \sum_{i=1}^{n} x_i f(x_i)
Variance:
\begin{aligned} \sigma^2 = V(X) &= E[(X - \mu)^2] \\ &= \sum_{i=1}^{n} (x_i - \mu)^2 f(x_i) \\ &= \sum_{i=1}^{n} x_i^2 f(x_i) - \mu^2 \end{aligned}
Standard Deviation: \sigma = \sqrt{V(X)}
Interpretation:
For the transmission error example:
Mean:
\begin{aligned} \mu &= E(X) = 0(0.6561) + 1(0.2916) + 2(0.0486) + 3(0.0036) + 4(0.0001) \\ &= 0 + 0.2916 + 0.0972 + 0.0108 + 0.0004 = 0.4 \end{aligned}
Variance:
First, calculate E(X^2):
E(X^2) = 0^2(0.6561) + 1^2(0.2916) + 2^2(0.0486) + 3^2(0.0036) + 4^2(0.0001) = 0.52
Then: V(X) = E(X^2) - [E(X)]^2 = 0.52 - (0.4)^2 = 0.52 - 0.16 = 0.36
Standard Deviation: \sigma = \sqrt{0.36} = 0.6
Error in eval(e, .data, pe): object 'e_x_squared' not found[1] "Mean and Variance for Transmission Errors:"Error: object 'discrete_stats' not foundError in eval(ei, .data, pe): object 'discrete_stats' not found[1] "Detailed Variance Calculation:"Error: object 'calc_table' not found# Example 3-14: Mean and Variance Calculation for Discrete Distribution
discrete_stats <- transmission_data %>%
fmutate(
x_times_p = errors * probability,
x_squared_times_p = errors^2 * probability
) %>%
fsummarise(
mean_x = fsum(x_times_p),
e_x_squared = fsum(x_squared_times_p),
variance_x = e_x_squared - mean_x^2,
sd_x = sqrt(variance_x)
)
print("Mean and Variance for Transmission Errors:")
print(discrete_stats)
# Create detailed calculation table
calc_table <- transmission_data %>%
fmutate(
x_times_p = errors * probability,
x_minus_mu = errors - discrete_stats$mean_x[1],
x_minus_mu_squared = x_minus_mu^2,
variance_term = x_minus_mu_squared * probability
)
print("Detailed Variance Calculation:")
print(calc_table)Two product designs are compared based on revenue potential:
Analysis:
Design A: E(X_A) = 3 million dollars, V(X_A) = 0
Design B: \begin{aligned} E(X_B) &= 7(0.3) + 2(0.7) = 2.1 + 1.4 = 3.5 \text{ million} \\ V(X_B) &= (7-3.5)^2(0.3) + (2-3.5)^2(0.7) \\ &= (3.5)^2(0.3) + (-1.5)^2(0.7) \\ &= 3.675 + 1.575 = 5.25 \text{ million}^2 \end{aligned}
Decision: Design B has higher expected revenue but much higher risk (variance).
A binomial experiment consists of n independent trials, each with:
Definition: The random variable X that counts the number of successes in n binomial trials has a binomial distribution with parameters n and p.
Probability Mass Function:
f(x) = P(X = x) = \binom{n}{x} p^x (1-p)^{n-x}
for x = 0, 1, 2, \ldots, n, where \binom{n}{x} = \frac{n!}{x!(n-x)!}.
Parameters:
Mean and Variance:
\begin{aligned} E(X) &= np \\ V(X) &= np(1-p) \end{aligned}
Notation: X \sim \text{Binomial}(n, p) or X \sim B(n, p)
Each water sample has a 10% chance of containing high levels of organic solids. Assume samples are independent. For the next 18 samples:
Problem 1: Probability that exactly 2 contain high solids?
Solution:
X \sim B(18, 0.1)
P(X = 2) = \binom{18}{2} (0.1)^2 (0.9)^{16} = 153 \times 0.01 \times 0.1853 = 0.284
Problem 2: Probability that at least 4 samples contain high solids?
Solution:
\begin{aligned} P(X \geq 4) &= 1 - P(X \leq 3) \\ &= 1 - \sum_{x=0}^{3} \binom{18}{x} (0.1)^x (0.9)^{18-x} \\ &= 1 - [0.150 + 0.300 + 0.284 + 0.168] = 0.098 \end{aligned}
Problem 3: Mean and variance?
Solution:
\begin{aligned} E(X) &= np = 18 \times 0.1 = 1.8 \\ V(X) &= np(1-p) = 18 \times 0.1 \times 0.9 = 1.62 \end{aligned}
[1] "Water Sample Binomial Analysis:" n p test_values prob_exactly_2 prob_at_least_4 prob_at_least_4_alt
<num> <num> <list> <num> <num> <num>
1: 18 0.1 2,4 0.2835121 0.09819684 0.09819684
prob_3_to_7 mean_binomial variance_binomial sd_binomial
<num> <num> <num> <num>
1: 0.2660305 1.8 1.62 1.272792[1] "Binomial Probability Table (first 9 values):" x probability cumulative
<int> <num> <num>
1: 0 0.1500946353 0.1500946
2: 1 0.3001892706 0.4502839
3: 2 0.2835120889 0.7337960
4: 3 0.1680071638 0.9018032
5: 4 0.0700029849 0.9718061
6: 5 0.0217787064 0.9935848
7: 6 0.0052430219 0.9988279
8: 7 0.0009986708 0.9998265
9: 8 0.0001525747 0.9999791# Example 3-16: Water Sample Analysis (Binomial)
water_sample_data <- data.table(
n = 18,
p = 0.1,
test_values = list(c(2, 4))
) %>%
fmutate(
# Problem 1: Exactly 2 samples
prob_exactly_2 = dbinom(2, size = n, prob = p),
# Problem 2: At least 4 samples
prob_at_least_4 = 1 - pbinom(3, size = n, prob = p),
prob_at_least_4_alt = pbinom(3, size = n, prob = p, lower.tail = FALSE),
# Problem 3: Between 3 and 7 (inclusive)
prob_3_to_7 = pbinom(7, size = n, prob = p) - pbinom(2, size = n, prob = p),
# Mean and variance
mean_binomial = n * p,
variance_binomial = n * p * (1 - p),
sd_binomial = sqrt(variance_binomial)
)
print("Water Sample Binomial Analysis:")
print(water_sample_data)
# Create binomial probability table
binomial_table <- data.table(
x = 0:8
) %>%
fmutate(
probability = dbinom(x, size = 18, prob = 0.1),
cumulative = pbinom(x, size = 18, prob = 0.1)
)
print("Binomial Probability Table (first 9 values):")
print(binomial_table)
fwrite(binomial_table, "data/binomial_probabilities.csv")A Poisson process models events occurring randomly over time or space. Examples include:
Definition: A Poisson process with rate \lambda has these properties:
Applications:
If events occur according to a Poisson process with rate \lambda, then the number of events X in a unit interval has a Poisson distribution.
Probability Mass Function:
f(x) = P(X = x) = \frac{\lambda^x e^{-\lambda}}{x!}
for x = 0, 1, 2, \ldots
Parameter:
Mean and Variance:
E(X) = V(X) = \lambda
Key Property: Mean equals variance (equidispersion)
Notation: X \sim \text{Poisson}(\lambda)
Flaws in copper wire follow a Poisson distribution with mean 2.3 flaws per millimeter.
Problem 1: Probability of exactly 2 flaws in 1 millimeter?
Solution: X \sim \text{Poisson}(2.3), P(X = 2) = \frac{e^{-2.3} \times 2.3^2}{2!} = \frac{e^{-2.3} \times 5.29}{2} = 0.265
Problem 2: Probability of 10 flaws in 5 millimeters?
Solution: For 5 mm, \lambda = 5 \times 2.3 = 11.5 P(X = 10) = \frac{e^{-11.5} \times 11.5^{10}}{10!} = 0.113
[1] "Poisson Distribution Analysis:" lambda_per_mm test_scenarios prob_2_in_1mm lambda_5mm prob_10_in_5mm
<num> <list> <num> <num> <num>
1: 2.3 1,5 0.2651846 11.5 0.1129351
mean_1mm variance_1mm mean_5mm variance_5mm
<num> <num> <num> <num>
1: 2.3 2.3 11.5 11.5[1] "Poisson Table for 1mm (lambda = 2.3):" x probability cumulative
<int> <num> <num>
1: 0 0.10025884 0.1002588
2: 1 0.23059534 0.3308542
3: 2 0.26518464 0.5960388
4: 3 0.20330823 0.7993471
5: 4 0.11690223 0.9162493
6: 5 0.05377503 0.9700243# Example 3-17: Copper Wire Flaws (Poisson)
poisson_data <- data.table(
lambda_per_mm = 2.3,
test_scenarios = list(c(1, 5)) # 1 mm and 5 mm
) %>%
fmutate(
# Problem 1: Exactly 2 flaws in 1 mm
prob_2_in_1mm = dpois(2, lambda = lambda_per_mm),
# Problem 2: 10 flaws in 5 mm
lambda_5mm = lambda_per_mm * 5,
prob_10_in_5mm = dpois(10, lambda = lambda_5mm),
# Additional calculations
mean_1mm = lambda_per_mm,
variance_1mm = lambda_per_mm,
mean_5mm = lambda_5mm,
variance_5mm = lambda_5mm
)
print("Poisson Distribution Analysis:")
print(poisson_data)
# Create Poisson probability tables
poisson_1mm <- data.table(x = 0:10) %>%
fmutate(
probability = dpois(x, lambda = 2.3),
cumulative = ppois(x, lambda = 2.3)
)
poisson_5mm <- data.table(x = 0:20) %>%
fmutate(
probability = dpois(x, lambda = 11.5),
cumulative = ppois(x, lambda = 11.5)
)
print("Poisson Table for 1mm (lambda = 2.3):")
print(head(poisson_1mm))
fwrite(poisson_1mm, "data/poisson_1mm.csv")
fwrite(poisson_5mm, "data/poisson_5mm.csv")The exponential distribution models the time between events in a Poisson process.
Probability Density Function:
f(x) = \lambda e^{-\lambda x}, \quad x \geq 0
Cumulative Distribution Function:
F(x) = 1 - e^{-\lambda x}, \quad x \geq 0
Parameter:
Mean and Variance:
\begin{aligned} E(X) &= \frac{1}{\lambda} \\ V(X) &= \frac{1}{\lambda^2} \end{aligned}
Memoryless Property: P(X > s + t | X > s) = P(X > t)
Applications: Reliability analysis, queuing systems, survival analysis
User log-ons to a system follow a Poisson process with mean 25 log-ons per hour.
Problem: Probability of no log-ons in a 6-minute interval?
Solution: Let X = time until first log-on (in hours). Then X \sim \text{Exponential}(25).
We want P(X > 0.1) since 6 minutes = 0.1 hour.
P(X > 0.1) = e^{-25 \times 0.1} = e^{-2.5} = 0.082
Alternative approach: Number of events in 0.1 hour ~ Poisson(2.5) P(\text{0 events}) = \frac{e^{-2.5} \times 2.5^0}{0!} = e^{-2.5} = 0.082
[1] "Network Log-on Analysis:" lambda_per_hour time_interval_minutes time_interval_hours prob_no_logons
<num> <num> <num> <num>
1: 25 6 0.1 0.082085
prob_no_logons_alt expected_events prob_0_events_poisson mean_time_between
<num> <num> <num> <num>
1: 0.082085 2.5 0.082085 0.04
mean_time_minutes
<num>
1: 2.4[1] "Exponential Probabilities for Various Times:" time_hours prob_exceed prob_within time_minutes
<num> <num> <num> <num>
1: 0.1 8.208500e-02 0.9179150 6
2: 0.2 6.737947e-03 0.9932621 12
3: 0.5 3.726653e-06 0.9999963 30
4: 1.0 1.388794e-11 1.0000000 60# Example 3-18: Computer Network Log-ons (Exponential)
network_data <- data.table(
lambda_per_hour = 25,
time_interval_minutes = 6
) %>%
fmutate(
time_interval_hours = time_interval_minutes / 60,
# Probability of no log-ons in 6 minutes
prob_no_logons = pexp(time_interval_hours, rate = lambda_per_hour, lower.tail = FALSE),
prob_no_logons_alt = exp(-lambda_per_hour * time_interval_hours),
# Verification using Poisson (events in interval)
expected_events = lambda_per_hour * time_interval_hours,
prob_0_events_poisson = dpois(0, lambda = expected_events),
# Mean time between log-ons
mean_time_between = 1 / lambda_per_hour,
# Convert to minutes
mean_time_minutes = mean_time_between * 60
)
print("Network Log-on Analysis:")
print(network_data)
# Additional exponential calculations
exp_calculations <- data.table(
time_hours = c(0.1, 0.2, 0.5, 1.0)
) %>%
fmutate(
prob_exceed = pexp(time_hours, rate = 25, lower.tail = FALSE),
prob_within = pexp(time_hours, rate = 25),
time_minutes = time_hours * 60
)
print("Exponential Probabilities for Various Times:")
print(exp_calculations)When n is large, the binomial distribution approaches normality due to the Central Limit Theorem.
Approximation: If X \sim B(n, p), then for large n: Z = \frac{X - np}{\sqrt{np(1-p)}} \approx N(0, 1)
Rule of Thumb: Use when np > 5 and n(1-p) > 5
Continuity Correction: Since we’re approximating discrete with continuous:
where Y \sim N(np, np(1-p)).
For n = 50 bits transmitted with error probability p = 0.1:
Exact calculation: P(X \leq 2) = \sum_{x=0}^{2} \binom{50}{x} (0.1)^x (0.9)^{50-x} = 0.1117
Normal approximation: X \sim B(50, 0.1) approximately N(5, 4.5)
With continuity correction:
\begin{aligned} P(X \leq 2) &\approx P\left(Z \leq \frac{2.5 - 5}{\sqrt{4.5}}\right) \\ &= P(Z \leq -1.18) = 0.119 \end{aligned}
Comparison: Exact = 0.1117, Approximation = 0.119 (very close!)
[1] "Normal Approximation to Binomial:" n p test_value np nq conditions_met exact_prob mu_approx
<num> <num> <num> <num> <num> <lgcl> <num> <num>
1: 50 0.1 2 5 45 FALSE 0.1117288 5
sigma_approx z_score_no_cc approx_no_cc z_score_with_cc approx_with_cc
<num> <num> <num> <num> <num>
1: 2.12132 -1.414214 0.0786496 -1.178511 0.1192964
error_no_cc error_with_cc
<num> <num>
1: 0.03307915 0.007567658[1] "Comparison Table (first 9 values):" x exact approx_cc error
<int> <num> <num> <num>
1: 0 0.005153775 0.01694743 0.011793652
2: 1 0.033785860 0.04948008 0.015694217
3: 2 0.111728756 0.11929641 0.007567658
4: 3 0.250293906 0.23975006 0.010543845
5: 4 0.431198407 0.40683186 0.024366549
6: 5 0.616123008 0.59316814 0.022954866
7: 6 0.770226842 0.76024994 0.009976903
8: 7 0.877854916 0.88070359 0.002848669
9: 8 0.942132794 0.95051992 0.008387129# Example 3-19: Normal Approximation to Binomial
normal_approx_data <- data.table(
n = 50,
p = 0.1,
test_value = 2
) %>%
fmutate(
# Check conditions for normal approximation
np = n * p,
nq = n * (1 - p),
conditions_met = (np > 5) & (nq > 5),
# Exact binomial calculation
exact_prob = pbinom(test_value, size = n, prob = p),
# Normal approximation parameters
mu_approx = np,
sigma_approx = sqrt(np * (1 - p)),
# Normal approximation without continuity correction
z_score_no_cc = (test_value - mu_approx) / sigma_approx,
approx_no_cc = pnorm(z_score_no_cc),
# Normal approximation with continuity correction
z_score_with_cc = (test_value + 0.5 - mu_approx) / sigma_approx,
approx_with_cc = pnorm(z_score_with_cc),
# Error calculations
error_no_cc = abs(exact_prob - approx_no_cc),
error_with_cc = abs(exact_prob - approx_with_cc)
)
print("Normal Approximation to Binomial:")
print(normal_approx_data)
# Compare multiple values
comparison_table <- data.table(x = 0:8) %>%
fmutate(
exact = pbinom(x, size = 50, prob = 0.1),
approx_cc = pnorm((x + 0.5 - 5) / sqrt(4.5)),
error = abs(exact - approx_cc)
)
print("Comparison Table (first 9 values):")
print(comparison_table)For large \lambda, the Poisson distribution approaches normality.
Approximation: If X \sim \text{Poisson}(\lambda), then for large \lambda: Z = \frac{X - \lambda}{\sqrt{\lambda}} \approx N(0, 1)
Rule of Thumb: Use when \lambda > 5
Continuity Correction: Apply the same corrections as for binomial approximation.
Contamination particles in water follow Poisson(1000). Find P(X < 950).
Exact: P(X < 950) = \sum_{x=0}^{949} \frac{e^{-1000} \times 1000^x}{x!} (computationally intensive)
Normal approximation: X \sim \text{Poisson}(1000) approximately N(1000, 1000)
With continuity correction:
\begin{aligned} P(X < 950) &= P(X \leq 949) \\ &\approx P\left(Z \leq \frac{949.5 - 1000}{\sqrt{1000}}\right) \\ &= P(Z \leq -1.60) = 0.055 \end{aligned}
[1] "Normal Approximation to Poisson:" lambda test_value condition_met exact_prob_approx mu_approx sigma_approx
<num> <num> <lgcl> <num> <num> <num>
1: 1000 950 TRUE 0.05420667 1000 31.62278
z_score normal_approx approx_error
<num> <num> <num>
1: -1.59695 0.0551384 0.0009317283[1] "Poisson Approximation Examples:" lambda test_x exact z_score approx error
<num> <num> <num> <num> <num> <num>
1: 10 15 0.9512596 1.7392527 0.9590048 0.007745243
2: 25 30 0.8633089 1.1000000 0.8643339 0.001025070
3: 50 55 0.7844704 0.7778175 0.7816617 0.002808718
4: 100 105 0.7128079 0.5500000 0.7088403 0.003967569# Example 3-20: Normal Approximation to Poisson
poisson_approx_data <- data.table(
lambda = 1000,
test_value = 950
) %>%
fmutate(
# Check condition for normal approximation
condition_met = lambda > 5,
# Exact Poisson (approximated due to computational limits)
exact_prob_approx = ppois(test_value - 1, lambda = lambda), # P(X < 950) = P(X <= 949)
# Normal approximation parameters
mu_approx = lambda,
sigma_approx = sqrt(lambda),
# Normal approximation with continuity correction
z_score = (test_value - 0.5 - mu_approx) / sigma_approx,
normal_approx = pnorm(z_score),
# Error (approximate since exact is computationally intensive)
approx_error = abs(exact_prob_approx - normal_approx)
)
print("Normal Approximation to Poisson:")
print(poisson_approx_data)
# Additional Poisson approximation examples
poisson_examples <- data.table(
lambda = c(10, 25, 50, 100),
test_x = c(15, 30, 55, 105)
) %>%
fmutate(
exact = ppois(test_x, lambda = lambda),
z_score = (test_x + 0.5 - lambda) / sqrt(lambda),
approx = pnorm(z_score),
error = abs(exact - approx)
)
print("Poisson Approximation Examples:")
print(poisson_examples)When multiple random variables are measured simultaneously, we need joint probability distributions.
Joint PDF (Continuous): For continuous random variables X and Y:
P(a < X < b, c < Y < d) = \int_a^b \int_c^d f(x,y) \, dy \, dx
Properties:
Joint PMF (Discrete): For discrete random variables: f(x,y) = P(X = x, Y = y)
Marginal Distributions:
Random variables X_1, X_2, \ldots, X_n are independent if:
P(X_1 \in E_1, X_2 \in E_2, \ldots, X_n \in E_n) = P(X_1 \in E_1) \times P(X_2 \in E_2) \times \cdots \times P(X_n \in E_n)
for any sets E_1, E_2, \ldots, E_n.
Equivalent Conditions:
Practical Meaning: Knowledge about one variable provides no information about others.
Shaft diameters are normally distributed with mean 0.2508 inch and standard deviation 0.0005 inch. Each shaft has probability 0.919 of meeting specifications.
Problem: If diameters are independent, what’s the probability that all 10 shafts meet specifications?
Solution: \begin{aligned} P(\text{all 10 meet specs}) &= P(X_1 \text{ meets}) \times P(X_2 \text{ meets}) \times \cdots \times P(X_{10} \text{ meets}) \\ &= (0.919)^{10} = 0.430 \end{aligned}
Interpretation: Only 43% chance that all 10 shafts will be acceptable.
A system operates if there’s a functional path from left to right. Components function independently with given probabilities.
Solution: For a series system: P(\text{system works}) = P(C_1 \text{ and } C_2) = P(C_1) \times P(C_2) = 0.9 \times 0.95 = 0.855
For parallel systems, calculate probability that at least one component works.
For independent random variables X_1, X_2, \ldots, X_n with means \mu_i and variances \sigma_i^2, consider the linear function:
Y = c_0 + c_1X_1 + c_2X_2 + \cdots + c_nX_n
Mean and Variance:
\begin{aligned} E(Y) &= c_0 + c_1\mu_1 + c_2\mu_2 + \cdots + c_n\mu_n \\ V(Y) &= c_1^2\sigma_1^2 + c_2^2\sigma_2^2 + \cdots + c_n^2\sigma_n^2 \end{aligned}
Key Properties:
Special Case: If X_1, X_2, \ldots, X_n are independent and normally distributed, then Y is also normally distributed.
A rectangular part has length X_1 and width X_2 that are independent and normally distributed: - X_1 \sim N(2, 0.1^2) cm - X_2 \sim N(5, 0.2^2) cm
The perimeter is Y = 2X_1 + 2X_2.
Solution:
\begin{aligned} E(Y) &= 2E(X_1) + 2E(X_2) = 2(2) + 2(5) = 14 \text{ cm} \\ V(Y) &= 2^2V(X_1) + 2^2V(X_2) = 4(0.1^2) + 4(0.2^2) \\ &= 4(0.01) + 4(0.04) = 0.04 + 0.16 = 0.20 \text{ cm}^2 \\ \sigma_Y &= \sqrt{0.20} = 0.447 \text{ cm} \end{aligned}
Since X_1 and X_2 are normal and independent, Y \sim N(14, 0.20).
Probability calculation: P(Y > 14.5) = P\left(Z > \frac{14.5 - 14}{0.447}\right) = P(Z > 1.12) = 0.131
When random variables are not independent, we need to account for their relationship.
Covariance: \text{Cov}(X, Y) = E[(X - \mu_X)(Y - \mu_Y)] = E(XY) - \mu_X\mu_Y
Correlation Coefficient: \rho_{XY} = \frac{\text{Cov}(X, Y)}{\sigma_X \sigma_Y}
Properties:
General Variance Formula: For Y = c_0 + c_1X_1 + c_2X_2 + \cdots + c_nX_n:
V(Y) = \sum_{i=1}^n c_i^2\sigma_i^2 + 2\sum_{i<j} c_ic_j\text{Cov}(X_i, X_j)
An investment portfolio consists of two stocks with returns X_1 and X_2:
Portfolio return: Y = 0.6X_1 + 0.4X_2
Solution:
\begin{aligned} E(Y) &= 0.6(0.12) + 0.4(0.08) = 0.072 + 0.032 = 0.104 \\ V(Y) &= (0.6)^2(0.04) + (0.4)^2(0.02) + 2(0.6)(0.4)(0.01) \\ &= 0.0144 + 0.0032 + 0.0048 = 0.0224 \end{aligned}
Interpretation: Positive covariance increases portfolio risk compared to independent case.
For nonlinear functions, we use Taylor series approximation around the means.
Single Variable: If Y = h(X) where X has mean \mu_X and variance \sigma_X^2:
\begin{aligned} E(Y) &\approx h(\mu_X) \\ V(Y) &\approx \left(\frac{dh}{dx}\bigg|_{x=\mu_X}\right)^2 \sigma_X^2 \end{aligned}
Multiple Variables: If Y = h(X_1, X_2, \ldots, X_n) with independent X_i:
\begin{aligned} E(Y) &\approx h(\mu_1, \mu_2, \ldots, \mu_n) \\ V(Y) &\approx \sum_{i=1}^n \left(\frac{\partial h}{\partial x_i}\bigg|_{\boldsymbol{\mu}}\right)^2 \sigma_i^2 \end{aligned}
Applications: Engineering tolerance analysis, measurement uncertainty
Power dissipated by a resistor: P = I^2R where I is current and R is resistance.
Given: I \sim N(20, 0.1^2) amperes, R = 80 ohms (constant)
Solution: h(I) = I^2R = 80I^2, so \frac{dh}{dI} = 160I
\begin{aligned} E(P) &\approx h(\mu_I) = 80(20)^2 = 32,000 \text{ watts} \\ V(P) &\approx \left(\frac{dh}{dI}\bigg|_{I=20}\right)^2 \sigma_I^2 \\ &= (160 \times 20)^2 (0.1)^2 = (3200)^2(0.01) = 102,400 \text{ watts}^2 \\ \sigma_P &= \sqrt{102,400} = 320 \text{ watts} \end{aligned}
[1] "Power Dissipation Analysis:" current_mean current_sd resistance derivative_at_mean power_mean_approx
<num> <num> <num> <num> <num>
1: 20 0.1 80 3200 32000
power_variance_approx power_sd_approx exact_mean exact_variance exact_sd
<num> <num> <num> <num> <num>
1: 102400 320 32000.8 102401.3 320.002[1] "Multiple Variable Function Analysis:" mu1 sigma1_sq mu2 sigma2_sq mu3 sigma3_sq dY_dX1 dY_dX2 dY_dX3
<num> <num> <num> <num> <num> <num> <num> <num> <num>
1: 10 1 5 0.25 2 0.04 2.5 5 -12.5
Y_mean_approx Y_variance_approx Y_sd_approx
<num> <num> <num>
1: 25 18.75 4.330127# Example 3-25: Electrical Power Dissipation (Propagation of Error)
power_analysis <- data.table(
current_mean = 20,
current_sd = 0.1,
resistance = 80
) %>%
fmutate(
# Function: P = I^2 * R
# Derivative: dP/dI = 2*I*R
derivative_at_mean = 2 * current_mean * resistance,
# Mean of P
power_mean_approx = current_mean^2 * resistance,
# Variance of P using propagation of error
power_variance_approx = (derivative_at_mean)^2 * current_sd^2,
power_sd_approx = sqrt(power_variance_approx),
# Exact calculation for comparison (since I is normal, I^2 follows chi-square scaled)
# For normal I, E[I^2] = mu^2 + sigma^2, Var[I^2] = 2*sigma^4 + 4*mu^2*sigma^2
exact_mean = (current_mean^2 + current_sd^2) * resistance,
exact_variance = (2 * current_sd^4 + 4 * current_mean^2 * current_sd^2) * resistance^2,
exact_sd = sqrt(exact_variance)
)
print("Power Dissipation Analysis:")
print(power_analysis)
# Multiple variable example: Y = X1*X2/X3
multi_var_data <- data.table(
mu1 = 10, sigma1_sq = 1,
mu2 = 5, sigma2_sq = 0.25,
mu3 = 2, sigma3_sq = 0.04
) %>%
fmutate(
# Function: Y = X1*X2/X3
# Partial derivatives at means
dY_dX1 = mu2 / mu3,
dY_dX2 = mu1 / mu3,
dY_dX3 = -(mu1 * mu2) / (mu3^2),
# Approximate mean
Y_mean_approx = (mu1 * mu2) / mu3,
# Approximate variance
Y_variance_approx = (dY_dX1)^2 * sigma1_sq + (dY_dX2)^2 * sigma2_sq + (dY_dX3)^2 * sigma3_sq,
Y_sd_approx = sqrt(Y_variance_approx)
)
print("Multiple Variable Function Analysis:")
print(multi_var_data)Consider Y = \frac{X_1X_2}{X_3} where X_1, X_2, X_3 are independent with:
Solution:
\frac{\partial h}{\partial x_1} = \frac{x_2}{x_3}, \quad \frac{\partial h}{\partial x_2} = \frac{x_1}{x_3}, \quad \frac{\partial h}{\partial x_3} = -\frac{x_1x_2}{x_3^2}
At (\mu_1, \mu_2, \mu_3) = (10, 5, 2):
\frac{\partial h}{\partial x_1}\bigg|_{\boldsymbol{\mu}} = \frac{5}{2} = 2.5, \quad \frac{\partial h}{\partial x_2}\bigg|_{\boldsymbol{\mu}} = \frac{10}{2} = 5, \quad \frac{\partial h}{\partial x_3}\bigg|_{\boldsymbol{\mu}} = -\frac{50}{4} = -12.5
\begin{aligned} E(Y) &\approx \frac{10 \times 5}{2} = 25 \\ V(Y) &\approx (2.5)^2(1) + (5)^2(0.25) + (-12.5)^2(0.04) \\ &= 6.25 + 6.25 + 6.25 = 18.75 \end{aligned}
Random Sample: Independent random variables X_1, X_2, \ldots, X_n with the same distribution are called a random sample of size n.
Statistic: A statistic is any function of the random variables in a random sample. Examples include:
Sampling Distribution: The probability distribution of a statistic is called its sampling distribution.
Key Properties of Sample Mean: If X_1, X_2, \ldots, X_n is a random sample from a population with mean \mu and variance \sigma^2:
\begin{aligned} E(\overline{X}) &= \mu \\ V(\overline{X}) &= \frac{\sigma^2}{n} \\ \sigma_{\overline{X}} &= \frac{\sigma}{\sqrt{n}} \end{aligned}
Central Limit Theorem: If X_1, X_2, \ldots, X_n is a random sample of size n from a population with mean \mu and variance \sigma^2, then as n \to \infty:
Z = \frac{\overline{X} - \mu}{\sigma/\sqrt{n}} \xrightarrow{d} N(0, 1)
Key Points:
Practical Importance: - Justifies normal approximations - Foundation for confidence intervals - Basis for hypothesis testing - Explains why sample means are approximately normal
An electronics company manufactures resistors with mean resistance 100 Ω and standard deviation 10 Ω. Find the probability that a random sample of n = 25 resistors has average resistance less than 95 Ω.
Solution: Let \overline{X} be the sample mean resistance.
By the Central Limit Theorem: \overline{X} \sim N\left(100, \frac{10^2}{25}\right) = N(100, 4)
\begin{aligned} P(\overline{X} < 95) &= P\left(\frac{\overline{X} - 100}{10/\sqrt{25}} < \frac{95 - 100}{2}\right) \\ &= P(Z < -2.5) = 0.0062 \end{aligned}
Interpretation: Only 0.62% chance that the sample mean is below 95 Ω.
[1] "Central Limit Theorem - Resistor Example:" population_mean population_sd sample_size test_value sample_mean_mean
<num> <num> <num> <num> <num>
1: 100 10 25 95 100
sample_mean_sd sample_mean_variance z_score probability percentile_5
<num> <num> <num> <num> <num>
1: 2 4 -2.5 0.006209665 96.71029
percentile_95
<num>
1: 103.2897[1] "CLT Demonstration - Effect of Sample Size:" sample_size sample_mean_sd prob_below_95 prob_above_105
<num> <num> <num> <num>
1: 1 10.000000 3.085375e-01 3.085375e-01
2: 4 5.000000 1.586553e-01 1.586553e-01
3: 9 3.333333 6.680720e-02 6.680720e-02
4: 16 2.500000 2.275013e-02 2.275013e-02
5: 25 2.000000 6.209665e-03 6.209665e-03
6: 36 1.666667 1.349898e-03 1.349898e-03
7: 49 1.428571 2.326291e-04 2.326291e-04
8: 64 1.250000 3.167124e-05 3.167124e-05# Example 3-27: Resistor Manufacturing (Central Limit Theorem)
resistor_clt <- data.table(
population_mean = 100,
population_sd = 10,
sample_size = 25,
test_value = 95
) %>%
fmutate(
# Sample mean distribution
sample_mean_mean = population_mean,
sample_mean_sd = population_sd / sqrt(sample_size),
sample_mean_variance = sample_mean_sd^2,
# Probability calculation
z_score = (test_value - sample_mean_mean) / sample_mean_sd,
probability = pnorm(z_score),
# Additional percentiles
percentile_5 = qnorm(0.05, mean = sample_mean_mean, sd = sample_mean_sd),
percentile_95 = qnorm(0.95, mean = sample_mean_mean, sd = sample_mean_sd)
)
print("Central Limit Theorem - Resistor Example:")
print(resistor_clt)
# CLT demonstration with different sample sizes
clt_demo <- data.table(
sample_size = c(1, 4, 9, 16, 25, 36, 49, 64)
) %>%
fmutate(
sample_mean_sd = 10 / sqrt(sample_size),
prob_below_95 = pnorm(95, mean = 100, sd = sample_mean_sd),
prob_above_105 = pnorm(105, mean = 100, sd = sample_mean_sd, lower.tail = FALSE)
)
print("CLT Demonstration - Effect of Sample Size:")
print(clt_demo)A production process has mean output 500 units/hour with standard deviation 50 units/hour. If we observe production for 36 hours, what’s the probability that average hourly production exceeds 510 units?
Solution: n = 36, \mu = 500, \sigma = 50
\begin{aligned} P(\overline{X} > 510) &= P\left(Z > \frac{510 - 500}{50/\sqrt{36}}\right) \\ &= P\left(Z > \frac{10}{50/6}\right) \\ &= P(Z > 1.2) = 1 - 0.8849 = 0.1151 \end{aligned}
Interpretation: About 11.5% chance of observing such high average production.
Continuous Random Variables:
Important Continuous Distributions:
Discrete Random Variables:
Important Discrete Distributions:
Multiple Random Variables:
Central Limit Theorem:
Applications:
[1] "Distribution Summary Statistics:" distribution sample_mean sample_sd sample_min sample_max sample_range
<char> <num> <num> <num> <num> <num>
1: Normal 50.90406 9.128159 26.90831124 71.87333 44.96502
2: Exponential 10.43276 9.353413 0.04367371 43.65911 43.61544
3: Weibull 88.00317 42.605619 5.82109597 189.09828 183.27719
4: Lognormal 60.70836 31.776841 10.42973266 186.75246 176.32273
cv
<num>
1: 0.1793208
2: 0.8965424
3: 0.4841373
4: 0.5234344[1] "Parameter Estimates:" distribution param1_name param1_estimate param2_name param2_estimate
<char> <char> <num> <char> <num>
1: Normal mean 50.9040591 sd 9.0824033
2: Exponential rate 0.0958519 <NA> NA
3: Weibull shape 2.1616282 scale 98.9731545
4: Lognormal meanlog 3.9722223 sdlog 0.5320079[1] "Goodness of Fit Tests:" distribution ks_statistic ks_p_value good_fit
<char> <num> <num> <lgcl>
1: Normal 0.05871913 0.8807748 TRUE
2: Exponential 0.07511794 0.6251819 TRUE
3: Weibull 0.09167937 0.3699697 TRUE
4: Lognormal 0.09690997 0.3046243 TRUE# Comprehensive Summary and Distribution Fitting Examples
library(fitdistrplus)
# Generate sample data for distribution fitting
set.seed(123)
# Example datasets
normal_sample <- rnorm(100, mean = 50, sd = 10)
exponential_sample <- rexp(100, rate = 0.1)
weibull_sample <- rweibull(100, shape = 2, scale = 100)
lognormal_sample <- rlnorm(100, meanlog = 4, sdlog = 0.5)
# Create comprehensive dataset
sample_data <- data.table(
normal_data = normal_sample,
exponential_data = exponential_sample,
weibull_data = weibull_sample,
lognormal_data = lognormal_sample
) %>%
fmutate(
observation = 1:fnobs(normal_data)
)
# Fit distributions using fitdistrplus
fit_normal <- fitdist(sample_data$normal_data, "norm")
fit_exponential <- fitdist(sample_data$exponential_data, "exp")
fit_weibull <- fitdist(sample_data$weibull_data, "weibull")
fit_lognormal <- fitdist(sample_data$lognormal_data, "lnorm")
# Summary statistics for each distribution
distribution_summary <- data.table(
distribution = c("Normal", "Exponential", "Weibull", "Lognormal"),
sample_mean = c(
fmean(sample_data$normal_data),
fmean(sample_data$exponential_data),
fmean(sample_data$weibull_data),
fmean(sample_data$lognormal_data)
),
sample_sd = c(
fsd(sample_data$normal_data),
fsd(sample_data$exponential_data),
fsd(sample_data$weibull_data),
fsd(sample_data$lognormal_data)
),
sample_min = c(
fmin(sample_data$normal_data),
fmin(sample_data$exponential_data),
fmin(sample_data$weibull_data),
fmin(sample_data$lognormal_data)
),
sample_max = c(
fmax(sample_data$normal_data),
fmax(sample_data$exponential_data),
fmax(sample_data$weibull_data),
fmax(sample_data$lognormal_data)
)
) %>%
fmutate(
sample_range = sample_max - sample_min,
cv = sample_sd / sample_mean
)
print("Distribution Summary Statistics:")
print(distribution_summary)
# Parameter estimates
param_estimates <- data.table(
distribution = c("Normal", "Exponential", "Weibull", "Lognormal"),
param1_name = c("mean", "rate", "shape", "meanlog"),
param1_estimate = c(
fit_normal$estimate[1],
fit_exponential$estimate[1],
fit_weibull$estimate[1],
fit_lognormal$estimate[1]
),
param2_name = c("sd", NA, "scale", "sdlog"),
param2_estimate = c(
fit_normal$estimate[2],
NA,
fit_weibull$estimate[2],
fit_lognormal$estimate[2]
)
)
print("Parameter Estimates:")
print(param_estimates)
# Goodness of fit tests
gof_tests <- data.table(
distribution = c("Normal", "Exponential", "Weibull", "Lognormal"),
ks_statistic = c(
ks.test(sample_data$normal_data, "pnorm",
mean = fit_normal$estimate[1], sd = fit_normal$estimate[2]
)$statistic,
ks.test(sample_data$exponential_data, "pexp",
rate = fit_exponential$estimate[1]
)$statistic,
ks.test(sample_data$weibull_data, "pweibull",
shape = fit_weibull$estimate[1], scale = fit_weibull$estimate[2]
)$statistic,
ks.test(sample_data$lognormal_data, "plnorm",
meanlog = fit_lognormal$estimate[1], sdlog = fit_lognormal$estimate[2]
)$statistic
),
ks_p_value = c(
ks.test(sample_data$normal_data, "pnorm",
mean = fit_normal$estimate[1], sd = fit_normal$estimate[2]
)$p.value,
ks.test(sample_data$exponential_data, "pexp",
rate = fit_exponential$estimate[1]
)$p.value,
ks.test(sample_data$weibull_data, "pweibull",
shape = fit_weibull$estimate[1], scale = fit_weibull$estimate[2]
)$p.value,
ks.test(sample_data$lognormal_data, "plnorm",
meanlog = fit_lognormal$estimate[1], sdlog = fit_lognormal$estimate[2]
)$p.value
)
) %>%
fmutate(
good_fit = ks_p_value > 0.05
)
print("Goodness of Fit Tests:")
print(gof_tests)
# Save all datasets
fwrite(sample_data, "data/distribution_samples.csv")
fwrite(distribution_summary, "data/distribution_summary.csv")
fwrite(param_estimates, "data/parameter_estimates.csv")
fwrite(gof_tests, "data/goodness_of_fit.csv"):::