Statistics is not only a perfect means for (empirical) research but also helpful in software engineering. This tutorial covers statistics for software engineers in a practical way
INTRODUCTION
Software architecture exists at the intersection of technical design, business requirements, and empirical validation. While architectural patterns and principles provide guidance, the effectiveness of architectural decisions must ultimately be validated through measurement and analysis. Statistics provides software architects with rigorous methods to evaluate system performance, understand user behavior, validate design choices, and make evidence-based decisions about scalability, reliability, and quality attributes.
This tutorial explores how statistical methods and empirical research can enhance software architecture practice. We will examine fundamental statistical concepts with deep explanations of the underlying mathematics, explore various testing methodologies with detailed algorithmic descriptions, and demonstrate practical applications through code examples. Every statistical term, technique, and algorithm will be thoroughly explained so that you understand not just how to use these tools, but why they work and what they truly mean.
PART ONE: FOUNDATIONAL STATISTICAL CONCEPTS
Understanding Populations and Samples
In software architecture, we frequently need to make inferences about system behavior based on limited observations. A population represents the complete set of all possible instances we wish to study. For example, the population might be all requests that will ever be processed by a service over its entire operational lifetime, all users who will interact with a system, or all possible execution paths through a codebase. The population is typically so large or infinite that we cannot measure it completely.
A sample consists of the actual measurements we collect from the population. When we run a load test and capture one thousand response times, those one thousand measurements constitute a sample. The fundamental challenge of statistics is using information from the sample to make valid inferences about the population.
The relationship between sample and population determines the validity of our conclusions. If our sample systematically differs from the population in important ways, we have sampling bias. For instance, if we only test our system during off-peak hours, our sample will underestimate the response times users actually experience. Understanding and controlling sampling bias is critical for architectural decision-making.
The Arithmetic Mean: Understanding Averages
The arithmetic mean, commonly called the average, is the most fundamental measure of central tendency. To calculate the mean, we sum all values in our dataset and divide by the count of values. Mathematically, if we have n observations denoted as x1, x2, x3, up to xn, the mean is calculated as the sum of all xi divided by n.
The mean represents the balance point of the data. If you imagine each data point as a weight on a number line, the mean is the point where the line would balance perfectly. This geometric interpretation helps us understand why the mean is sensitive to extreme values: a single very large value shifts the balance point substantially.
For software architects analyzing response times, the mean provides a useful summary but can be misleading when the distribution is skewed. If most requests complete in 50 milliseconds but occasional requests take 5000 milliseconds, the mean will be inflated by those slow requests and may not represent typical user experience.
Let us examine an implementation that calculates the mean and explores its properties:
# Arithmetic mean calculation and analysis
def calculate_mean(data):
"""
Calculate the arithmetic mean (average) of a dataset.
The mean is the sum of all values divided by the count.
Formula: mean = (x1 + x2 + ... + xn) / n
Parameters:
data: List of numerical values
Returns:
Float representing the mean, or None if data is empty
"""
if not data or len(data) == 0:
return None
total = sum(data)
count = len(data)
mean = total / count
return mean
Here is the weighted mean implementation:
def calculate_weighted_mean(values, weights):
"""
Calculate weighted mean where each value has an associated weight.
Weighted mean gives more influence to values with higher weights.
Formula: weighted_mean = sum(xi * wi) / sum(wi)
This is useful when different observations have different importance
or represent different sample sizes.
Parameters:
values: List of numerical values
weights: List of weights corresponding to each value
Returns:
Float representing the weighted mean
"""
if len(values) != len(weights):
raise ValueError("Values and weights must have same length")
if not values:
return None
weighted_sum = sum(v * w for v, w in zip(values, weights))
total_weight = sum(weights)
if total_weight == 0:
raise ValueError("Total weight cannot be zero")
return weighted_sum / total_weight
Now let us see an example usage:
# Example: Response time analysis
response_times = [45, 52, 48, 51, 49, 53, 47, 50, 48, 51]
mean_response = calculate_mean(response_times)
print(f"Mean response time: {mean_response:.2f} ms")
# Demonstrate sensitivity to outliers
response_times_with_outlier = response_times + [5000]
mean_with_outlier = calculate_mean(response_times_with_outlier)
print(f"Mean with outlier: {mean_with_outlier:.2f} ms")
print(f"Impact of single outlier: {mean_with_outlier - mean_response:.2f} ms increase")
# Weighted mean example: combining metrics from different servers
# Each server processed different numbers of requests
server_means = [45.2, 48.7, 46.3] # Mean response time per server
request_counts = [1000, 1500, 800] # Number of requests per server
overall_mean = calculate_weighted_mean(server_means, request_counts)
print(f"\nWeighted mean across servers: {overall_mean:.2f} ms")
This implementation demonstrates both simple and weighted means. The weighted mean is particularly important in distributed systems where you need to combine metrics from components handling different volumes of traffic. The code shows how a single outlier can dramatically affect the mean, illustrating why software architects must understand the distribution of their data, not just its average.
The Median: The Middle Value
The median represents the middle value when all observations are arranged in order from smallest to largest. Unlike the mean, which uses all values in its calculation, the median depends only on the position of values. This positional nature makes the median resistant to extreme values, which is why it often better represents typical behavior in skewed distributions.
To calculate the median, we first sort the data. If we have an odd number of observations, the median is the value in the middle position. If we have an even number of observations, the median is conventionally defined as the average of the two middle values. For example, with five values arranged as 10, 20, 30, 40, 50, the median is 30. With six values arranged as 10, 20, 30, 40, 50, 60, the median is 35, which is the average of 30 and 40.
The median divides the dataset into two equal halves: half the observations are below the median and half are above it. This property makes the median the 50th percentile, a concept we will explore in detail later. For software architects analyzing response times, the median often provides a better sense of typical user experience than the mean, especially when occasional very slow requests exist.
Here is an implementation that calculates the median and explains the algorithm:
# Median calculation with detailed explanation
def calculate_median(data):
"""
Calculate the median (middle value) of a dataset.
Algorithm:
1. Sort the data in ascending order
2. If odd number of values: return the middle value
3. If even number of values: return average of two middle values
The median is the 50th percentile, meaning 50% of values
are below it and 50% are above it.
Parameters:
data: List of numerical values
Returns:
Float representing the median
"""
if not data or len(data) == 0:
return None
# Step 1: Sort the data
# Sorting is O(n log n) complexity
sorted_data = sorted(data)
n = len(sorted_data)
# Step 2: Find the middle position(s)
if n % 2 == 1:
# Odd number of values: one middle position
middle_index = n // 2
median = sorted_data[middle_index]
else:
# Even number of values: two middle positions
upper_middle_index = n // 2
lower_middle_index = upper_middle_index - 1
# Average the two middle values
median = (sorted_data[lower_middle_index] +
sorted_data[upper_middle_index]) / 2
return median
Here is a function to compare mean and median:
def compare_mean_and_median(data):
"""
Compare mean and median to understand data distribution.
Relationship between mean and median reveals skewness:
- Mean > Median: Right-skewed (long tail of high values)
- Mean < Median: Left-skewed (long tail of low values)
- Mean approximately equal to Median: Symmetric distribution
Parameters:
data: List of numerical values
Returns:
Dictionary with mean, median, and interpretation
"""
mean = calculate_mean(data)
median = calculate_median(data)
difference = mean - median
if abs(difference) < 0.01 * mean:
skewness = "approximately symmetric"
elif difference > 0:
skewness = "right-skewed (mean pulled up by high values)"
else:
skewness = "left-skewed (mean pulled down by low values)"
return {
'mean': mean,
'median': median,
'difference': difference,
'skewness': skewness
}
Example usage:
# Example: Analyzing database query times
query_times = [12, 15, 13, 18, 14, 16, 15, 17, 14, 15,
13, 16, 14, 15, 180] # Note the outlier
analysis = compare_mean_and_median(query_times)
print("Database Query Time Analysis")
print("=" * 50)
print(f"Mean: {analysis['mean']:.2f} ms")
print(f"Median: {analysis['median']:.2f} ms")
print(f"Difference: {analysis['difference']:.2f} ms")
print(f"Distribution: {analysis['skewness']}")
print()
# Demonstrate robustness of median
normal_queries = [q for q in query_times if q < 100]
print(f"Without outlier - Mean: {calculate_mean(normal_queries):.2f} ms")
print(f"Without outlier - Median: {calculate_median(normal_queries):.2f} ms")
This implementation shows how the median resists the influence of outliers. The comparison function reveals the shape of the distribution by examining the relationship between mean and median. For software architects, this distinction is crucial: a service with median response time of 50 milliseconds but mean of 200 milliseconds has a serious tail latency problem that the median alone would not reveal.
Percentiles and Quantiles: Understanding Distribution Position
Percentiles generalize the concept of the median to arbitrary positions in the distribution. The pth percentile is the value below which p percent of the observations fall. The median is the 50th percentile because 50 percent of observations fall below it. The 95th percentile, commonly used in service level objectives, is the value below which 95 percent of observations fall.
Percentiles are also called quantiles when expressed as proportions rather than percentages. The 0.95 quantile is equivalent to the 95th percentile. Quartiles divide the data into four equal parts: the first quartile is the 25th percentile, the second quartile is the 50th percentile or median, and the third quartile is the 75th percentile.
Calculating percentiles requires careful handling of edge cases and interpolation. When the desired percentile does not correspond exactly to a data point, we must interpolate between adjacent values. Several interpolation methods exist, with linear interpolation being most common. The algorithm first calculates the position in the sorted data corresponding to the desired percentile, then interpolates if that position falls between two data points.
Here is a comprehensive implementation of percentile calculation:
# Percentile calculation with multiple interpolation methods
import math
def calculate_percentile(data, p, method='linear'):
"""
Calculate the pth percentile of a dataset.
Percentiles indicate the value below which a given percentage
of observations fall. For example, the 95th percentile is the
value below which 95% of the data lies.
Algorithm (linear interpolation method):
1. Sort the data
2. Calculate position: pos = (p/100) * (n-1)
3. If pos is integer: return data[pos]
4. If pos is fractional: interpolate between floor(pos) and ceil(pos)
Parameters:
data: List of numerical values
p: Percentile to calculate (0-100)
method: Interpolation method ('linear', 'lower', 'higher', 'nearest')
Returns:
Float representing the percentile value
"""
if not data or len(data) == 0:
return None
if p < 0 or p > 100:
raise ValueError("Percentile must be between 0 and 100")
sorted_data = sorted(data)
n = len(sorted_data)
# Handle edge cases
if p == 0:
return sorted_data[0]
if p == 100:
return sorted_data[-1]
# Calculate position in sorted array
# Using (n-1) makes 0th percentile = min and 100th percentile = max
position = (p / 100.0) * (n - 1)
if method == 'linear':
# Linear interpolation between adjacent values
lower_index = int(math.floor(position))
upper_index = int(math.ceil(position))
if lower_index == upper_index:
return sorted_data[lower_index]
# Interpolation fraction
fraction = position - lower_index
# Linear interpolation formula
lower_value = sorted_data[lower_index]
upper_value = sorted_data[upper_index]
result = lower_value + fraction * (upper_value - lower_value)
return result
elif method == 'lower':
# Always use the lower adjacent value
return sorted_data[int(math.floor(position))]
elif method == 'higher':
# Always use the higher adjacent value
return sorted_data[int(math.ceil(position))]
elif method == 'nearest':
# Use the nearest value
return sorted_data[int(round(position))]
else:
raise ValueError(f"Unknown method: {method}")
Here is the quartiles calculation:
def calculate_quartiles(data):
"""
Calculate the three quartiles that divide data into four equal parts.
Q1 (First Quartile): 25th percentile
Q2 (Second Quartile): 50th percentile (median)
Q3 (Third Quartile): 75th percentile
The interquartile range (IQR = Q3 - Q1) measures the spread
of the middle 50% of the data.
Parameters:
data: List of numerical values
Returns:
Dictionary with Q1, Q2, Q3, and IQR
"""
q1 = calculate_percentile(data, 25)
q2 = calculate_percentile(data, 50)
q3 = calculate_percentile(data, 75)
iqr = q3 - q1
return {
'Q1': q1,
'Q2': q2,
'Q3': q3,
'IQR': iqr
}
SLA percentiles function:
def calculate_sla_percentiles(data):
"""
Calculate common SLA percentiles for performance analysis.
Service Level Agreements often specify percentile-based targets:
- P50 (median): Typical user experience
- P95: Experience of 95% of users
- P99: Experience of 99% of users (tail latency)
- P99.9: Extreme tail latency
Parameters:
data: List of numerical values (e.g., response times)
Returns:
Dictionary with common percentiles
"""
return {
'P50': calculate_percentile(data, 50),
'P90': calculate_percentile(data, 90),
'P95': calculate_percentile(data, 95),
'P99': calculate_percentile(data, 99),
'P99.9': calculate_percentile(data, 99.9)
}
Example usage:
# Example: API response time percentile analysis
response_times = [
45, 48, 52, 47, 50, 49, 51, 48, 46, 53,
47, 50, 49, 52, 48, 51, 47, 49, 50, 48,
53, 48, 50, 52, 49, 51, 47, 54, 48, 50,
55, 58, 62, 68, 75, 82, 95, 110, 125, 180
]
print("API Response Time Percentile Analysis")
print("=" * 60)
# Calculate SLA percentiles
sla_metrics = calculate_sla_percentiles(response_times)
for percentile, value in sla_metrics.items():
print(f"{percentile}: {value:.2f} ms")
print()
# Calculate quartiles
quartiles = calculate_quartiles(response_times)
print("Quartile Analysis:")
print(f"Q1 (25th percentile): {quartiles['Q1']:.2f} ms")
print(f"Q2 (50th percentile/median): {quartiles['Q2']:.2f} ms")
print(f"Q3 (75th percentile): {quartiles['Q3']:.2f} ms")
print(f"IQR (Interquartile Range): {quartiles['IQR']:.2f} ms")
print()
print("Interpretation:")
print(f"- 50% of requests complete in {sla_metrics['P50']:.2f} ms or less")
print(f"- 95% of requests complete in {sla_metrics['P95']:.2f} ms or less")
print(f"- Only 1% of requests exceed {sla_metrics['P99']:.2f} ms")
This implementation provides multiple interpolation methods for percentile calculation. The linear interpolation method is most common and provides smooth percentile estimates. The code demonstrates how percentiles reveal the full distribution shape, not just central tendency. For software architects, percentiles are essential for understanding tail latency and setting realistic service level objectives.
Variance: Measuring Spread Around the Mean
Variance quantifies how spread out the data is around the mean. It measures the average squared deviation from the mean, providing a single number that captures the variability in the dataset. Understanding variance is fundamental to statistics because it appears in nearly every advanced statistical method.
To calculate variance, we first find the mean of the dataset. Then for each observation, we calculate its deviation from the mean by subtracting the mean from the observation. We square each deviation to ensure that positive and negative deviations do not cancel out, and to give greater weight to larger deviations. Finally, we average these squared deviations.
The formula for population variance uses n in the denominator, where n is the total count of observations. However, when calculating variance from a sample to estimate population variance, we use n minus one in the denominator. This adjustment, called Bessel's correction, compensates for the fact that sample variance tends to underestimate population variance. The reason is subtle: when we calculate sample variance, we use the sample mean rather than the true population mean, and the sample mean is calculated from the same data, introducing a dependency that reduces the apparent variability.
Here is a detailed implementation explaining variance calculation:
# Variance calculation with detailed mathematical explanation
def calculate_variance(data, sample=True):
"""
Calculate the variance of a dataset.
Variance measures the average squared deviation from the mean.
It quantifies how spread out the data is.
Mathematical formula:
Population variance: sigma squared = sum((xi - mu) squared) / n
Sample variance: s squared = sum((xi - x_bar) squared) / (n - 1)
Where:
- xi represents each individual observation
- mu is the population mean
- x_bar is the sample mean
- n is the number of observations
- sum means "sum of"
Why divide by (n-1) for sample variance?
This is Bessel's correction. When estimating population variance
from a sample, using n would systematically underestimate the
true variance. Using (n-1) provides an unbiased estimate.
The intuition: we "lose" one degree of freedom because we must
first estimate the mean from the data. Once we know n-1 values
and the mean, the nth value is determined.
Parameters:
data: List of numerical values
sample: If True, use sample variance (n-1); if False, use population variance (n)
Returns:
Float representing the variance
"""
if not data or len(data) == 0:
return None
n = len(data)
# Step 1: Calculate the mean
mean = sum(data) / n
# Step 2: Calculate squared deviations from mean
squared_deviations = []
for value in data:
deviation = value - mean
squared_deviation = deviation ** 2
squared_deviations.append(squared_deviation)
# Step 3: Sum the squared deviations
sum_squared_deviations = sum(squared_deviations)
# Step 4: Divide by n or (n-1)
if sample and n > 1:
# Sample variance: divide by (n-1)
variance = sum_squared_deviations / (n - 1)
else:
# Population variance: divide by n
variance = sum_squared_deviations / n
return variance
Alternative computational formula:
def calculate_variance_alternative(data, sample=True):
"""
Alternative variance calculation using computational formula.
This formula is mathematically equivalent but more efficient
for computation and less prone to rounding errors.
Computational formula:
s squared = [sum(xi squared) - (sum(xi)) squared / n] / (n-1)
This formula requires only one pass through the data if we
accumulate both sum(xi) and sum(xi squared) simultaneously.
Parameters:
data: List of numerical values
sample: If True, use sample variance
Returns:
Float representing the variance
"""
if not data or len(data) == 0:
return None
n = len(data)
# Calculate sum and sum of squares in one pass
sum_x = 0
sum_x_squared = 0
for value in data:
sum_x += value
sum_x_squared += value ** 2
# Apply computational formula
numerator = sum_x_squared - (sum_x ** 2) / n
denominator = (n - 1) if sample and n > 1 else n
variance = numerator / denominator
return variance
Demonstrating variance properties:
def demonstrate_variance_properties(data):
"""
Demonstrate key properties of variance.
Properties:
1. Variance is always non-negative
2. Variance is zero only when all values are identical
3. Variance has squared units (e.g., ms squared)
4. Adding a constant to all values doesn't change variance
5. Multiplying all values by constant c multiplies variance by c squared
Parameters:
data: List of numerical values
Returns:
Dictionary demonstrating properties
"""
original_variance = calculate_variance(data)
# Property 4: Adding constant doesn't change variance
shifted_data = [x + 100 for x in data]
shifted_variance = calculate_variance(shifted_data)
# Property 5: Scaling by c multiplies variance by c squared
scaled_data = [x * 2 for x in data]
scaled_variance = calculate_variance(scaled_data)
return {
'original_variance': original_variance,
'shifted_variance': shifted_variance,
'variance_unchanged_by_shift': abs(original_variance - shifted_variance) < 0.0001,
'scaled_variance': scaled_variance,
'scaling_factor_squared': (scaled_variance / original_variance) if original_variance > 0 else None
}
Example usage:
# Example: Response time variability analysis
response_times = [45, 52, 48, 51, 49, 53, 47, 50, 48, 51]
variance = calculate_variance(response_times, sample=True)
variance_alt = calculate_variance_alternative(response_times, sample=True)
print("Response Time Variability Analysis")
print("=" * 60)
print(f"Sample variance: {variance:.2f} ms squared")
print(f"Variance (computational formula): {variance_alt:.2f} ms squared")
print(f"Formulas agree: {abs(variance - variance_alt) < 0.0001}")
print()
# Demonstrate properties
properties = demonstrate_variance_properties(response_times)
print("Variance Properties:")
print(f"Original variance: {properties['original_variance']:.2f} ms squared")
print(f"Variance after adding 100: {properties['shifted_variance']:.2f} ms squared")
print(f"Unchanged by shift: {properties['variance_unchanged_by_shift']}")
print(f"Variance after multiplying by 2: {properties['scaled_variance']:.2f} ms squared")
print(f"Scaling factor: {properties['scaling_factor_squared']:.2f} (should be 4)")
This implementation explains both the definitional and computational formulas for variance. The code demonstrates that variance has squared units, which makes direct interpretation difficult. This limitation motivates the standard deviation, which we will explore next.
Standard Deviation: Variance in Original Units
Standard deviation is simply the square root of variance. This transformation returns the measure of spread to the original units of measurement. If we are measuring response times in milliseconds and calculate a variance of 25 square milliseconds, the standard deviation is 5 milliseconds. The standard deviation of 5 milliseconds is directly interpretable: it represents a typical deviation from the mean.
Standard deviation is denoted by the Greek letter sigma when referring to population standard deviation, and by the letter s when referring to sample standard deviation. The calculation is straightforward: compute the variance, then take its square root. All the properties of variance apply to standard deviation, with appropriate adjustments for the square root transformation.
For normally distributed data, which we will explore in depth later, standard deviation has a special interpretation. Approximately 68 percent of observations fall within one standard deviation of the mean, 95 percent within two standard deviations, and 99.7 percent within three standard deviations. This empirical rule, also called the 68-95-99.7 rule, provides intuition for interpreting standard deviation.
Here is an implementation with detailed explanation:
# Standard deviation calculation and interpretation
import math
def calculate_standard_deviation(data, sample=True):
"""
Calculate the standard deviation of a dataset.
Standard deviation is the square root of variance.
It measures spread in the same units as the original data.
Mathematical formula:
Population std dev: sigma = square root of [sum((xi - mu) squared) / n]
Sample std dev: s = square root of [sum((xi - x_bar) squared) / (n - 1)]
Standard deviation is preferred over variance for interpretation
because it has the same units as the data. For example, if
measuring response times in milliseconds, standard deviation
is also in milliseconds, while variance would be in ms squared.
Parameters:
data: List of numerical values
sample: If True, use sample std dev; if False, use population std dev
Returns:
Float representing the standard deviation
"""
variance = calculate_variance(data, sample)
if variance is None:
return None
# Standard deviation is the square root of variance
std_dev = math.sqrt(variance)
return std_dev
Coefficient of variation:
def calculate_coefficient_of_variation(data):
"""
Calculate the coefficient of variation (CV).
CV is the ratio of standard deviation to mean, expressed as a percentage.
Formula: CV = (s / x_bar) times 100%
The coefficient of variation is dimensionless, making it useful
for comparing variability across datasets with different units
or different scales. For example, you can compare the relative
variability of response times (measured in ms) to throughput
(measured in requests/sec).
Interpretation:
- CV < 10%: Low variability
- CV 10-20%: Moderate variability
- CV > 20%: High variability
Parameters:
data: List of numerical values
Returns:
Float representing CV as a percentage
"""
mean = calculate_mean(data)
std_dev = calculate_standard_deviation(data)
if mean == 0:
return float('inf')
cv = (std_dev / mean) * 100
return cv
Empirical rule analysis:
def empirical_rule_analysis(data):
"""
Apply the empirical rule (68-95-99.7 rule) to analyze data.
The empirical rule states that for normally distributed data:
- Approximately 68% of values fall within 1 std dev of mean
- Approximately 95% of values fall within 2 std dev of mean
- Approximately 99.7% of values fall within 3 std dev of mean
This function calculates these ranges and determines what
percentage of actual data falls within them.
Parameters:
data: List of numerical values
Returns:
Dictionary with empirical rule ranges and actual percentages
"""
mean = calculate_mean(data)
std_dev = calculate_standard_deviation(data)
n = len(data)
# Calculate ranges
one_sigma_lower = mean - std_dev
one_sigma_upper = mean + std_dev
two_sigma_lower = mean - 2 * std_dev
two_sigma_upper = mean + 2 * std_dev
three_sigma_lower = mean - 3 * std_dev
three_sigma_upper = mean + 3 * std_dev
# Count values within each range
within_one_sigma = sum(1 for x in data if one_sigma_lower <= x <= one_sigma_upper)
within_two_sigma = sum(1 for x in data if two_sigma_lower <= x <= two_sigma_upper)
within_three_sigma = sum(1 for x in data if three_sigma_lower <= x <= three_sigma_upper)
return {
'mean': mean,
'std_dev': std_dev,
'one_sigma_range': (one_sigma_lower, one_sigma_upper),
'one_sigma_percentage': (within_one_sigma / n) * 100,
'two_sigma_range': (two_sigma_lower, two_sigma_upper),
'two_sigma_percentage': (within_two_sigma / n) * 100,
'three_sigma_range': (three_sigma_lower, three_sigma_upper),
'three_sigma_percentage': (within_three_sigma / n) * 100
}
Example usage:
# Example: Analyzing service response time consistency
service_a_times = [45, 48, 52, 47, 50, 49, 51, 48, 46, 53]
service_b_times = [50, 85, 30, 65, 40, 75, 55, 45, 70, 35]
print("Service Response Time Consistency Analysis")
print("=" * 60)
# Service A analysis
mean_a = calculate_mean(service_a_times)
std_a = calculate_standard_deviation(service_a_times)
cv_a = calculate_coefficient_of_variation(service_a_times)
print("Service A:")
print(f" Mean: {mean_a:.2f} ms")
print(f" Standard Deviation: {std_a:.2f} ms")
print(f" Coefficient of Variation: {cv_a:.2f}%")
print()
# Service B analysis
mean_b = calculate_mean(service_b_times)
std_b = calculate_standard_deviation(service_b_times)
cv_b = calculate_coefficient_of_variation(service_b_times)
print("Service B:")
print(f" Mean: {mean_b:.2f} ms")
print(f" Standard Deviation: {std_b:.2f} ms")
print(f" Coefficient of Variation: {cv_b:.2f}%")
print()
print("Comparison:")
print(f"Service A has {cv_a:.1f}% relative variability")
print(f"Service B has {cv_b:.1f}% relative variability")
print(f"Service B is {cv_b/cv_a:.1f}x more variable than Service A")
print()
# Empirical rule analysis for Service A
empirical = empirical_rule_analysis(service_a_times)
print("Empirical Rule Analysis for Service A:")
print(f"1 std dev range: ({empirical['one_sigma_range'][0]:.2f}, {empirical['one_sigma_range'][1]:.2f}) ms")
print(f" Expected: ~68%, Actual: {empirical['one_sigma_percentage']:.1f}%")
print(f"2 std dev range: ({empirical['two_sigma_range'][0]:.2f}, {empirical['two_sigma_range'][1]:.2f}) ms")
print(f" Expected: ~95%, Actual: {empirical['two_sigma_percentage']:.1f}%")
This implementation demonstrates how standard deviation provides interpretable measures of variability. The coefficient of variation enables comparison of relative variability across different metrics or systems. The empirical rule analysis shows how standard deviation relates to the proportion of data within specific ranges, providing intuition for what a given standard deviation value means in practice.
PART TWO: PROBABILITY DISTRIBUTIONS
What is a Probability Distribution?
A probability distribution describes how probability is distributed across possible values of a variable. For software architects, probability distributions model the behavior of system metrics such as response times, error rates, or request arrival patterns. Understanding distributions enables us to calculate probabilities, set thresholds, and make predictions about system behavior.
Probability distributions come in two flavors: discrete and continuous. Discrete distributions apply to variables that can only take specific values, such as the number of errors in a time period or the count of concurrent users. Continuous distributions apply to variables that can take any value within a range, such as response time or CPU utilization percentage.
Every probability distribution is characterized by parameters that determine its shape and location. For example, the normal distribution has two parameters: mean and standard deviation. The exponential distribution has one parameter: the rate. Understanding these parameters and how they affect the distribution is essential for applying distributions to real-world problems.
The Normal Distribution: The Bell Curve
The normal distribution, also called the Gaussian distribution, is the most important probability distribution in statistics. It appears throughout nature and engineering due to the central limit theorem, which states that the sum or average of many independent random variables tends toward a normal distribution regardless of the underlying distributions.
The normal distribution is characterized by its symmetric bell shape. It has two parameters: mu which is the mean and determines the center of the distribution, and sigma which is the standard deviation and determines the spread. The probability density function, which gives the relative likelihood of different values, is highest at the mean and decreases symmetrically on both sides.
The mathematical formula for the normal probability density function is complex but important to understand. For a value x, the probability density is calculated as one divided by sigma times the square root of two pi, multiplied by e raised to the power of negative one half times the quantity x minus mu divided by sigma, all squared. This formula ensures that the total area under the curve equals one, as required for any probability distribution.
Let‘s create an implementation of the normal distribution:
# Normal distribution implementation with detailed mathematics
import math
class NormalDistribution:
"""
Implementation of the normal (Gaussian) probability distribution.
The normal distribution is characterized by two parameters:
- mu: The mean, which determines the center
- sigma: The standard deviation, which determines the spread
The probability density function (PDF) is:
f(x) = (1 / (sigma times square root of (2 times pi))) times exp(-((x-mu) squared) / (2 times sigma squared))
Where:
- pi approximately equals 3.14159
- e (Euler's number) approximately equals 2.71828
- exp(y) means e raised to the power y
- square root means square root
Key properties:
- Symmetric around the mean
- Mean = Median = Mode
- 68% of values within 1 std dev of mean
- 95% of values within 2 std dev of mean
- 99.7% of values within 3 std dev of mean
"""
def __init__(self, mean, std_dev):
"""
Initialize a normal distribution.
Parameters:
mean: The mean (mu) of the distribution
std_dev: The standard deviation (sigma) of the distribution
"""
if std_dev <= 0:
raise ValueError("Standard deviation must be positive")
self.mean = mean
self.std_dev = std_dev
self.variance = std_dev ** 2
def pdf(self, x):
"""
Calculate the probability density function at point x.
The PDF gives the relative likelihood of different values.
Note: For continuous distributions, PDF values can exceed 1.
The PDF must integrate to 1 over all possible values.
Algorithm:
1. Calculate the coefficient: 1 / (sigma times square root of (2 times pi))
2. Calculate the exponent: -((x-mu) squared) / (2 times sigma squared)
3. Multiply coefficient by e raised to exponent
Parameters:
x: The point at which to evaluate the PDF
Returns:
Float representing the probability density
"""
# Coefficient: 1 / (sigma times square root of (2 times pi))
coefficient = 1.0 / (self.std_dev * math.sqrt(2 * math.pi))
# Exponent: -((x-mu) squared) / (2 times sigma squared)
deviation = x - self.mean
exponent = -(deviation ** 2) / (2 * self.variance)
# PDF value
pdf_value = coefficient * math.exp(exponent)
return pdf_value
def cdf(self, x):
"""
Calculate the cumulative distribution function at point x.
The CDF gives the probability that a random variable is
less than or equal to x: P(X <= x)
For the normal distribution, the CDF is:
Phi(x) = (1/2) times [1 + erf((x-mu) / (sigma times square root of 2))]
Where erf is the error function, a special function that
represents the integral of the standard normal PDF.
Algorithm:
1. Standardize: z = (x - mu) / sigma
2. Apply error function: erf(z / square root of 2)
3. Transform: (1 + erf(z / square root of 2)) / 2
Parameters:
x: The point at which to evaluate the CDF
Returns:
Float representing the cumulative probability
"""
# Standardize to z-score
z = (x - self.mean) / self.std_dev
# Use error function to calculate CDF
# erf is available in Python's math module
cdf_value = 0.5 * (1.0 + math.erf(z / math.sqrt(2)))
return cdf_value
def inverse_cdf(self, probability):
"""
Calculate the inverse CDF (quantile function).
Given a probability p, find the value x such that P(X <= x) = p.
This is useful for calculating percentiles.
The inverse CDF has no closed-form solution for the normal
distribution, so we use a numerical approximation called the
Beasley-Springer-Moro algorithm.
Parameters:
probability: Probability value between 0 and 1
Returns:
Float representing the value at that probability
"""
if probability <= 0 or probability >= 1:
raise ValueError("Probability must be between 0 and 1")
# Use approximation for standard normal inverse CDF
z = self._inverse_standard_normal_cdf(probability)
# Transform from standard normal to this distribution
x = self.mean + z * self.std_dev
return x
def _inverse_standard_normal_cdf(self, p):
"""
Approximate the inverse CDF of standard normal distribution.
This uses the Beasley-Springer-Moro algorithm, which provides
accurate approximations for probabilities between 0 and 1.
The algorithm uses different polynomial approximations for
the central region (absolute value of (p - 0.5) < 0.42) and tail regions.
Parameters:
p: Probability (0 < p < 1)
Returns:
Float representing the z-score
"""
# Coefficients for central region
a = [2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637]
b = [-8.47351093090, 23.08336743743, -21.06224101826, 3.13082909833]
# Coefficients for tail region
c = [0.3374754822726147, 0.9761690190917186, 0.1607979714918209,
0.0276438810333863, 0.0038405729373609, 0.0003951896511919,
0.0000321767881768, 0.0000002888167364, 0.0000003960315187]
y = p - 0.5
if abs(y) < 0.42:
# Central region: use rational approximation
r = y * y
numerator = (((a[3] * r + a[2]) * r + a[1]) * r + a[0]) * y
denominator = (((b[3] * r + b[2]) * r + b[1]) * r + b[0]) * r + 1.0
z = numerator / denominator
else:
# Tail region: use logarithmic approximation
if y > 0:
r = math.log(-math.log(1.0 - p))
else:
r = math.log(-math.log(p))
z = c[0]
for i in range(1, len(c)):
z += c[i] * (r ** i)
if y < 0:
z = -z
return z
def probability_between(self, lower, upper):
"""
Calculate probability that value falls between lower and upper.
P(lower <= X <= upper) = CDF(upper) - CDF(lower)
Parameters:
lower: Lower bound
upper: Upper bound
Returns:
Float representing the probability
"""
return self.cdf(upper) - self.cdf(lower)
def z_score(self, x):
"""
Calculate the z-score (standardized value) for x.
The z-score represents how many standard deviations
x is away from the mean.
Formula: z = (x - mu) / sigma
Interpretation:
- z = 0: x equals the mean
- z = 1: x is 1 std dev above mean
- z = -1: x is 1 std dev below mean
Parameters:
x: Value to standardize
Returns:
Float representing the z-score
"""
return (x - self.mean) / self.std_dev
Example usage:
# Example: Modeling database query execution times
# Assume query times follow normal distribution with mean 120ms, std dev 15ms
query_dist = NormalDistribution(mean=120, std_dev=15)
print("Database Query Time Analysis (Normal Distribution)")
print("=" * 70)
print(f"Distribution: Normal(mu={query_dist.mean}, sigma={query_dist.std_dev})")
print()
# Calculate probabilities
prob_under_100 = query_dist.cdf(100)
prob_under_150 = query_dist.cdf(150)
prob_between = query_dist.probability_between(100, 150)
print("Probability Calculations:")
print(f"P(query time <= 100ms) = {prob_under_100:.4f} ({prob_under_100*100:.2f}%)")
print(f"P(query time <= 150ms) = {prob_under_150:.4f} ({prob_under_150*100:.2f}%)")
print(f"P(100ms <= query time <= 150ms) = {prob_between:.4f} ({prob_between*100:.2f}%)")
print()
# Calculate percentiles for SLA definition
print("Percentile Analysis for SLA:")
for p in [50, 90, 95, 99]:
value = query_dist.inverse_cdf(p / 100.0)
print(f" {p}th percentile: {value:.2f} ms")
print()
# Z-score analysis
test_values = [90, 120, 150, 180]
print("Z-Score Analysis:")
for value in test_values:
z = query_dist.z_score(value)
print(f" {value} ms: z = {z:.2f} ({abs(z):.2f} std devs from mean)")
This comprehensive implementation explains the mathematics behind the normal distribution. The code demonstrates probability calculations, percentile determination, and z-score standardization. For software architects, the normal distribution models many aggregate metrics and provides a foundation for statistical inference and hypothesis testing.
The Exponential Distribution: Modeling Time Between Events
The exponential distribution models the time between events in a Poisson process, where events occur continuously and independently at a constant average rate. In software architecture, the exponential distribution commonly models time between failures, time between request arrivals, or service time duration.
The exponential distribution has one parameter: lambda, which represents the rate of occurrence. If lambda equals 0.5 failures per hour, then on average we expect one failure every two hours. The mean of the exponential distribution equals one divided by lambda, so higher rates correspond to shorter average times.
The probability density function for the exponential distribution is lambda times e raised to the power of negative lambda times x, where x is the time value. This function is highest at x equals zero and decreases exponentially as x increases, reflecting the fact that short intervals are most likely while long intervals become progressively less probable.
The exponential distribution has a unique property called memorylessness. This means that the probability of an event occurring in the next time interval is independent of how much time has already elapsed. Mathematically, the probability that we wait more than s plus t time units, given that we have already waited s time units, equals the probability that we wait more than t time units from the start. This property makes the exponential distribution appropriate for modeling random failures but inappropriate for modeling wear-out failures or scheduled events.
See the detailed implementation:
# Exponential distribution implementation with mathematical details
import math
import random
class ExponentialDistribution:
"""
Implementation of the exponential probability distribution.
The exponential distribution models the time between events
occurring at a constant average rate in a Poisson process.
Parameter:
- lambda: The rate parameter (events per unit time)
The probability density function (PDF) is:
f(x) = lambda times exp(-lambda times x) for x >= 0
The cumulative distribution function (CDF) is:
F(x) = 1 - exp(-lambda times x) for x >= 0
Key properties:
- Mean = 1/lambda
- Variance = 1/(lambda squared)
- Standard deviation = 1/lambda
- Memoryless: P(X > s+t | X > s) = P(X > t)
The memoryless property means the distribution "forgets" the past.
If failures follow an exponential distribution, a component that
has run for 10 hours is no more or less likely to fail in the
next hour than a brand new component.
"""
def __init__(self, rate=None, mean=None):
"""
Initialize an exponential distribution.
Specify either the rate (lambda) or the mean (1/lambda).
Parameters:
rate: Event rate (lambda parameter, events per unit time)
mean: Mean time between events (1/lambda)
"""
if rate is None and mean is None:
raise ValueError("Must specify either rate or mean")
if rate is not None and mean is not None:
raise ValueError("Specify only one of rate or mean, not both")
if rate is not None:
if rate <= 0:
raise ValueError("Rate must be positive")
self.rate = rate
self.mean = 1.0 / rate
else:
if mean <= 0:
raise ValueError("Mean must be positive")
self.mean = mean
self.rate = 1.0 / mean
# Variance and standard deviation
self.variance = 1.0 / (self.rate ** 2)
self.std_dev = 1.0 / self.rate
def pdf(self, x):
"""
Calculate the probability density function at time x.
PDF formula: f(x) = lambda times exp(-lambda times x)
The PDF gives the relative likelihood of different wait times.
For the exponential distribution, the PDF is highest at x=0
and decreases exponentially.
Algorithm:
1. Check that x >= 0 (negative times are impossible)
2. Calculate lambda times exp(-lambda times x)
Parameters:
x: Time value (must be non-negative)
Returns:
Float representing the probability density
"""
if x < 0:
return 0.0
# f(x) = lambda times exp(-lambda times x)
pdf_value = self.rate * math.exp(-self.rate * x)
return pdf_value
def cdf(self, x):
"""
Calculate the cumulative distribution function at time x.
CDF formula: F(x) = 1 - exp(-lambda times x)
The CDF gives the probability that an event occurs by time x.
This represents P(X <= x), where X is the waiting time.
Derivation:
The CDF is the integral of the PDF from 0 to x:
F(x) = integral from 0 to x of lambda times exp(-lambda times t) dt
= [-exp(-lambda times t)] evaluated from 0 to x
= -exp(-lambda times x) - (-exp(0))
= 1 - exp(-lambda times x)
Parameters:
x: Time value
Returns:
Float representing cumulative probability
"""
if x < 0:
return 0.0
# F(x) = 1 - exp(-lambda times x)
cdf_value = 1.0 - math.exp(-self.rate * x)
return cdf_value
def survival_function(self, t):
"""
Calculate the survival function (reliability function) at time t.
Survival function: S(t) = P(X > t) = 1 - F(t) = exp(-lambda times t)
This gives the probability that no event occurs before time t.
In reliability engineering, this represents the probability
that a component survives beyond time t.
The survival function is also called the complementary CDF.
Parameters:
t: Time value
Returns:
Float representing survival probability
"""
if t < 0:
return 1.0
# S(t) = exp(-lambda times t)
survival = math.exp(-self.rate * t)
return survival
def hazard_rate(self):
"""
Calculate the hazard rate (instantaneous failure rate).
The hazard rate is defined as:
h(t) = f(t) / S(t) = [lambda times exp(-lambda times t)] / [exp(-lambda times t)] = lambda
For the exponential distribution, the hazard rate is constant
and equals lambda. This constant hazard rate is equivalent to the
memoryless property.
In reliability engineering, constant hazard rate characterizes
the "useful life" period where failures occur randomly, not
due to wear-out or infant mortality.
Returns:
Float representing the hazard rate (equals lambda)
"""
return self.rate
def inverse_cdf(self, probability):
"""
Calculate the inverse CDF (quantile function).
Given probability p, find time x such that P(X <= x) = p.
Derivation:
Start with: p = 1 - exp(-lambda times x)
Rearrange: exp(-lambda times x) = 1 - p
Take ln: -lambda times x = ln(1 - p)
Solve for x: x = -ln(1 - p) / lambda
This closed-form solution makes the exponential distribution
particularly tractable for analytical calculations.
Parameters:
probability: Probability value between 0 and 1
Returns:
Float representing the time value
"""
if probability <= 0 or probability >= 1:
raise ValueError("Probability must be between 0 and 1")
# x = -ln(1 - p) / lambda
x = -math.log(1.0 - probability) / self.rate
return x
def percentile(self, p):
"""
Calculate the pth percentile.
Convenience method that converts percentile (0-100) to
probability (0-1) and calls inverse_cdf.
Parameters:
p: Percentile (0-100)
Returns:
Float representing the time value at that percentile
"""
if p < 0 or p > 100:
raise ValueError("Percentile must be between 0 and 100")
return self.inverse_cdf(p / 100.0)
def generate_sample(self, size=1):
"""
Generate random samples from this exponential distribution.
Uses the inverse transform method:
1. Generate uniform random number U between 0 and 1
2. Apply inverse CDF: X = -ln(1-U) / lambda
This method works because if U is uniform on (0,1), then
F inverse of (U) has distribution F.
Parameters:
size: Number of samples to generate
Returns:
List of random samples
"""
samples = []
for _ in range(size):
# Python's random.expovariate(lambda) generates exponential samples
sample = random.expovariate(self.rate)
samples.append(sample)
return samples
Example usage:
# Example: Analyzing service failure patterns
# Assume failures occur at average rate of 1 per 720 hours (30 days)
failure_dist = ExponentialDistribution(mean=720)
print("Service Reliability Analysis (Exponential Distribution)")
print("=" * 70)
print(f"Mean Time Between Failures (MTBF): {failure_dist.mean:.2f} hours")
print(f"Failure rate (lambda): {failure_dist.rate:.6f} per hour")
print(f"Standard deviation: {failure_dist.std_dev:.2f} hours")
print()
# Reliability calculations
print("Reliability Analysis:")
time_periods = [24, 168, 720, 8760] # 1 day, 1 week, 1 month, 1 year
for hours in time_periods:
reliability = failure_dist.survival_function(hours)
failure_prob = 1 - reliability
days = hours / 24
print(f" {days:5.0f} days: P(survive) = {reliability:.4f}, P(fail) = {failure_prob:.4f}")
print()
# Percentile analysis
print("Time to Failure Percentiles:")
for p in [10, 25, 50, 75, 90]:
time_hours = failure_dist.percentile(p)
time_days = time_hours / 24
print(f" {p:2d}th percentile: {time_days:6.1f} days ({time_hours:7.1f} hours)")
print()
# Demonstrate memoryless property
print("Memoryless Property Demonstration:")
print("If system has already run 30 days without failure:")
conditional_prob = failure_dist.survival_function(30 * 24)
unconditional_prob = failure_dist.survival_function(30 * 24)
print(f" P(survive next 30 days | survived 30 days) = {conditional_prob:.4f}")
print(f" P(survive 30 days from new) = {unconditional_prob:.4f}")
print(f" These are equal due to memoryless property")
This implementation thoroughly explains the exponential distribution's mathematics and unique memoryless property. The code demonstrates reliability analysis, percentile calculations, and the relationship between rate and mean. For software architects, the exponential distribution is essential for modeling random failures and analyzing system reliability.
The Poisson Distribution: Counting Events
The Poisson distribution models the number of events occurring in a fixed interval of time or space when events occur independently at a constant average rate. While the exponential distribution models time between events, the Poisson distribution models the count of events. These two distributions are intimately related: if inter-arrival times follow an exponential distribution with rate lambda, then the count of arrivals in a unit time interval follows a Poisson distribution with parameter lambda.
The Poisson distribution has one parameter: lambda, which represents the expected number of events in the interval. If we expect an average of 5 requests per second, lambda equals 5. The distribution gives probabilities for observing 0 events, 1 event, 2 events, and so on.
The probability mass function for the Poisson distribution gives the probability of observing exactly k events. The formula is lambda raised to the power k, times e raised to the power negative lambda, divided by k factorial. The factorial term k factorial equals k times k minus one times k minus two, continuing down to 1. For example, 5 factorial equals 5 times 4 times 3 times 2 times 1, which equals 120.
The Poisson distribution has an interesting property: its mean equals its variance, both equal to lambda. This relationship provides a diagnostic test: if we observe count data where the variance substantially exceeds the mean, the Poisson distribution may not be appropriate, suggesting overdispersion.
The following implementation demonstrates this:
# Poisson distribution implementation with detailed mathematics
import math
class PoissonDistribution:
"""
Implementation of the Poisson probability distribution.
The Poisson distribution models the number of events occurring
in a fixed interval when events occur independently at a
constant average rate.
Parameter:
- lambda: Expected number of events in the interval
The probability mass function (PMF) is:
P(X = k) = (lambda to the power k times e to the power negative lambda) / k!
Where:
- k is the number of events (0, 1, 2, ...)
- e approximately equals 2.71828 (Euler's number)
- k! is k factorial (k times (k-1) times (k-2) times ... times 1)
Key properties:
- Mean = lambda
- Variance = lambda
- Standard deviation = square root of lambda
- Mean equals variance (equidispersion)
Applications in software architecture:
- Request arrival rates
- Error counts in time windows
- Number of concurrent users
- Event occurrences in monitoring
"""
def __init__(self, lambda_param):
"""
Initialize a Poisson distribution.
Parameters:
lambda_param: Expected number of events (lambda, must be positive)
"""
if lambda_param <= 0:
raise ValueError("Lambda parameter must be positive")
self.lambda_param = lambda_param
self.mean = lambda_param
self.variance = lambda_param
self.std_dev = math.sqrt(lambda_param)
def pmf(self, k):
"""
Calculate probability mass function for exactly k events.
PMF formula: P(X = k) = (lambda to the power k times e to the power negative lambda) / k!
Algorithm:
1. Calculate lambda to the power k
2. Calculate e to the power negative lambda
3. Calculate k! (k factorial)
4. Multiply and divide: (lambda to the power k times e to the power negative lambda) / k!
Note: For large k, k! becomes very large and can cause
numerical overflow. We use logarithms to avoid this:
ln(PMF) = k times ln(lambda) - lambda - ln(k!)
PMF = exp(ln(PMF))
Parameters:
k: Number of events (must be non-negative integer)
Returns:
Float representing probability of exactly k events
"""
if k < 0:
raise ValueError("k must be non-negative")
if not isinstance(k, int):
raise ValueError("k must be an integer")
# Use logarithms to avoid numerical overflow
# ln(PMF) = k times ln(lambda) - lambda - ln(k!)
log_pmf = k * math.log(self.lambda_param) - self.lambda_param - math.lgamma(k + 1)
# Note: math.lgamma(k+1) = ln(k!)
# We use lgamma instead of factorial to handle large k
pmf_value = math.exp(log_pmf)
return pmf_value
def cdf(self, k):
"""
Calculate cumulative distribution function for at most k events.
CDF formula: P(X <= k) = sum from i=0 to k of P(X = i)
The CDF sums the probabilities of 0, 1, 2, ..., k events.
For the Poisson distribution, the CDF can also be expressed
using the incomplete gamma function, but direct summation
is straightforward for reasonable values of k.
Parameters:
k: Maximum number of events
Returns:
Float representing probability of k or fewer events
"""
if k < 0:
return 0.0
# Sum PMF values from 0 to k
cumulative = 0.0
for i in range(int(k) + 1):
cumulative += self.pmf(i)
return cumulative
def probability_at_least(self, k):
"""
Calculate probability of at least k events.
P(X >= k) = 1 - P(X <= k-1) = 1 - CDF(k-1)
This is the complement of the CDF at k-1.
Parameters:
k: Minimum number of events
Returns:
Float representing probability of k or more events
"""
if k <= 0:
return 1.0
return 1.0 - self.cdf(k - 1)
def probability_between(self, lower, upper):
"""
Calculate probability of event count between lower and upper.
P(lower <= X <= upper) = CDF(upper) - CDF(lower-1)
Parameters:
lower: Minimum number of events (inclusive)
upper: Maximum number of events (inclusive)
Returns:
Float representing the probability
"""
if lower > upper:
return 0.0
if lower == 0:
return self.cdf(upper)
return self.cdf(upper) - self.cdf(lower - 1)
def mode(self):
"""
Calculate the mode (most likely value).
For Poisson distribution:
- If lambda is an integer: modes are lambda and lambda-1 (bimodal)
- If lambda is not an integer: mode is floor(lambda)
Returns:
Integer or tuple representing the mode(s)
"""
if self.lambda_param == int(self.lambda_param):
# Lambda is an integer: bimodal
return (int(self.lambda_param), int(self.lambda_param) - 1)
else:
# Lambda is not an integer: single mode
return int(math.floor(self.lambda_param))
Example usage:
# Example: Analyzing API request arrival patterns
# Assume average of 25 requests per second
request_dist = PoissonDistribution(lambda_param=25)
print("API Request Arrival Analysis (Poisson Distribution)")
print("=" * 70)
print(f"Expected requests per second (lambda): {request_dist.lambda_param:.2f}")
print(f"Mean: {request_dist.mean:.2f}")
print(f"Variance: {request_dist.variance:.2f}")
print(f"Standard deviation: {request_dist.std_dev:.2f}")
print(f"Mode: {request_dist.mode()}")
print()
# Probability of specific counts
print("Probability of Specific Request Counts:")
for k in [15, 20, 25, 30, 35]:
prob = request_dist.pmf(k)
print(f" P(exactly {k} requests) = {prob:.4f} ({prob*100:.2f}%)")
print()
# Capacity planning: probability of exceeding thresholds
print("Capacity Planning - Probability of Exceeding Thresholds:")
thresholds = [30, 35, 40, 50]
for threshold in thresholds:
prob_exceed = request_dist.probability_at_least(threshold)
print(f" P({threshold}+ requests) = {prob_exceed:.4f} ({prob_exceed*100:.2f}%)")
print()
# Normal operating range
print("Normal Operating Range Analysis:")
prob_normal = request_dist.probability_between(20, 30)
print(f" P(20-30 requests) = {prob_normal:.4f} ({prob_normal*100:.2f}%)")
prob_low = request_dist.cdf(15)
print(f" P(<=15 requests, unusually low) = {prob_low:.4f} ({prob_low*100:.2f}%)")
prob_high = request_dist.probability_at_least(40)
print(f" P(>=40 requests, unusually high) = {prob_high:.4f} ({prob_high*100:.2f}%)")
print()
# Demonstrate mean = variance property
print("Poisson Property: Mean = Variance")
print(f" Mean: {request_dist.mean:.4f}")
print(f" Variance: {request_dist.variance:.4f}")
print(f" Ratio: {request_dist.variance / request_dist.mean:.4f} (should be 1.0)")
This implementation thoroughly explains the Poisson distribution's mathematics and properties. The code demonstrates probability calculations for capacity planning, anomaly detection, and understanding normal operating ranges. For software architects, the Poisson distribution is essential for modeling request arrivals, error counts, and other discrete event occurrences.
PART THREE: HYPOTHESIS TESTING
What is Hypothesis Testing?
Hypothesis testing provides a rigorous framework for making decisions based on empirical evidence while quantifying uncertainty. Software architects use hypothesis testing to evaluate whether observed differences between systems, configurations, or implementations represent genuine effects or merely random variation. The methodology protects against drawing conclusions from noise while providing a measure of confidence in our decisions.
The fundamental logic of hypothesis testing involves proof by contradiction. We begin by assuming that there is no effect, no difference, or no relationship. This assumption is called the null hypothesis, denoted H0. We then examine our data to see if it is compatible with this null hypothesis. If the data would be very unlikely under the null hypothesis, we reject the null hypothesis and conclude that an effect exists. If the data is reasonably compatible with the null hypothesis, we fail to reject it, meaning we do not have sufficient evidence to conclude an effect exists.
The alternative hypothesis, denoted H1 or Ha, represents the effect we seek to detect. For example, if we are comparing response times between two database configurations, the null hypothesis might state that the configurations have identical mean response times, while the alternative hypothesis states that they differ.
The p-value quantifies how compatible our data is with the null hypothesis. Specifically, the p-value is the probability of observing data at least as extreme as what we actually observed, assuming the null hypothesis is true. A small p-value indicates that our observed data would be unlikely if the null hypothesis were true, providing evidence against the null hypothesis.
We compare the p-value to a predetermined significance level, denoted alpha, which represents our threshold for rejecting the null hypothesis. Common choices are alpha equals 0.05 or 0.01. If the p-value is less than alpha, we reject the null hypothesis and conclude that the effect is statistically significant. If the p-value exceeds alpha, we fail to reject the null hypothesis.
This decision framework controls the Type I error rate, which is the probability of incorrectly rejecting a true null hypothesis. By setting alpha to 0.05, we ensure that if the null hypothesis is actually true, we will incorrectly reject it only 5 percent of the time in repeated experiments. Type II error represents the complementary mistake of failing to reject a false null hypothesis. The power of a test, defined as one minus the Type II error rate, represents the probability of correctly detecting an effect when it truly exists.
The T-Test: Comparing Means
The t-test enables comparison of means between groups when the underlying distributions are approximately normal. The test accounts for both the difference in means and the variability within each group, recognizing that the same mean difference is more convincing when variability is low than when it is high.
The test statistic, denoted t, measures how many standard errors the observed difference in means is from zero. The formula is the difference in sample means divided by the standard error of the difference. The standard error accounts for variability in both samples and their sizes.
For independent samples, where observations in one group are unrelated to observations in the other group, we use the independent samples t-test. This applies when comparing response times from two different server configurations, where each configuration processes different requests.
For paired samples, where observations come in matched pairs, we use the paired samples t-test. This applies when comparing before-and-after measurements on the same requests, or when comparing two methods applied to the same test cases. The paired test is more powerful because it eliminates variability between pairs.
The degrees of freedom, denoted df, depend on sample sizes and affect the critical values used to assess significance. For independent samples with equal variances, df equals n1 plus n2 minus 2. For unequal variances, we use the Welch-Satterthwaite approximation for degrees of freedom.
Below you’ll find the implementation:
# T-test implementation with detailed statistical theory
import math
class TTest:
"""
Implementation of t-tests for comparing means between two groups.
The t-test evaluates whether the difference between two sample
means is statistically significant, accounting for variability
and sample size.
Test statistic:
t = (x_bar_1 - x_bar_2) / SE
Where:
- x_bar_1, x_bar_2 are sample means
- SE is the standard error of the difference
The t-statistic follows a t-distribution under the null hypothesis.
The t-distribution is similar to the normal distribution but has
heavier tails, accounting for uncertainty in estimating variance
from small samples.
Types of t-tests:
1. Independent samples: Two separate groups
2. Paired samples: Matched observations
3. One-sample: Compare sample mean to known value
"""
def __init__(self, sample1, sample2=None):
"""
Initialize with sample data.
Parameters:
sample1: List of numerical values from first group
sample2: List of numerical values from second group (None for one-sample test)
"""
if len(sample1) < 2:
raise ValueError("Sample 1 must contain at least 2 observations")
self.sample1 = sample1
self.sample2 = sample2
self.n1 = len(sample1)
if sample2 is not None:
if len(sample2) < 2:
raise ValueError("Sample 2 must contain at least 2 observations")
self.n2 = len(sample2)
def _mean(self, data):
"""Calculate mean of a dataset."""
return sum(data) / len(data)
def _variance(self, data):
"""
Calculate sample variance.
Formula: s squared = sum((xi - x_bar) squared) / (n - 1)
Uses (n-1) for Bessel's correction to provide unbiased
estimate of population variance.
"""
mean = self._mean(data)
squared_devs = [(x - mean) ** 2 for x in data]
return sum(squared_devs) / (len(data) - 1)
def _std_dev(self, data):
"""Calculate sample standard deviation."""
return math.sqrt(self._variance(data))
def independent_samples_test(self, equal_variance=True):
"""
Perform independent samples t-test.
Tests whether two independent groups have different means.
Null hypothesis: mu_1 = mu_2 (population means are equal)
Alternative hypothesis: mu_1 not equal to mu_2 (population means differ)
Two variants:
1. Equal variance (Student's t-test): Assumes sigma_1 squared = sigma_2 squared
2. Unequal variance (Welch's t-test): Does not assume equal variance
Algorithm for equal variance:
1. Calculate pooled variance: sp squared = [(n_1-1)s_1 squared + (n_2-1)s_2 squared] / (n_1+n_2-2)
2. Calculate standard error: SE = square root of [sp squared times (1/n_1 + 1/n_2)]
3. Calculate t-statistic: t = (x_bar_1 - x_bar_2) / SE
4. Degrees of freedom: df = n_1 + n_2 - 2
Algorithm for unequal variance (Welch's t-test):
1. Calculate standard error: SE = square root of (s_1 squared/n_1 + s_2 squared/n_2)
2. Calculate t-statistic: t = (x_bar_1 - x_bar_2) / SE
3. Calculate Welch-Satterthwaite df:
df = (s_1 squared/n_1 + s_2 squared/n_2) squared / [(s_1 squared/n_1) squared/(n_1-1) + (s_2 squared/n_2) squared/(n_2-1)]
Parameters:
equal_variance: If True, use pooled variance; if False, use Welch's test
Returns:
Dictionary containing test results
"""
if self.sample2 is None:
raise ValueError("Independent samples test requires two samples")
# Calculate sample statistics
mean1 = self._mean(self.sample1)
mean2 = self._mean(self.sample2)
var1 = self._variance(self.sample1)
var2 = self._variance(self.sample2)
mean_diff = mean1 - mean2
if equal_variance:
# Pooled variance t-test (Student's t-test)
# Pooled variance: weighted average of the two variances
pooled_var = ((self.n1 - 1) * var1 + (self.n2 - 1) * var2) / (self.n1 + self.n2 - 2)
# Standard error using pooled variance
standard_error = math.sqrt(pooled_var * (1/self.n1 + 1/self.n2))
# Degrees of freedom
df = self.n1 + self.n2 - 2
test_type = "Student's t-test (equal variances)"
else:
# Welch's t-test (unequal variances)
# Standard error without pooling
standard_error = math.sqrt(var1/self.n1 + var2/self.n2)
# Welch-Satterthwaite degrees of freedom
numerator = (var1/self.n1 + var2/self.n2) ** 2
denominator = (var1/self.n1)**2 / (self.n1-1) + (var2/self.n2)**2 / (self.n2-1)
df = numerator / denominator
test_type = "Welch's t-test (unequal variances)"
# Calculate t-statistic
if standard_error == 0:
t_statistic = 0
else:
t_statistic = mean_diff / standard_error
# Calculate approximate p-value
p_value = self._approximate_t_test_p_value(t_statistic, df)
# Calculate confidence interval for difference
t_critical = self._t_critical_value(df, 0.05)
margin_of_error = t_critical * standard_error
ci_lower = mean_diff - margin_of_error
ci_upper = mean_diff + margin_of_error
return {
'test_type': test_type,
'mean1': mean1,
'mean2': mean2,
'mean_difference': mean_diff,
'std_error': standard_error,
't_statistic': t_statistic,
'degrees_of_freedom': df,
'p_value': p_value,
'significant_at_0.05': p_value < 0.05,
'significant_at_0.01': p_value < 0.01,
'confidence_interval_95': (ci_lower, ci_upper)
}
def paired_samples_test(self):
"""
Perform paired samples t-test.
Tests whether the mean difference between paired observations
is significantly different from zero.
Null hypothesis: mu_d = 0 (mean difference is zero)
Alternative hypothesis: mu_d not equal to 0 (mean difference is not zero)
Algorithm:
1. Calculate differences: di = x1i - x2i for each pair
2. Calculate mean difference: d_bar = sum(di) / n
3. Calculate standard deviation of differences: sd
4. Calculate standard error: SE = sd / square root of n
5. Calculate t-statistic: t = d_bar / SE
6. Degrees of freedom: df = n - 1
The paired test is more powerful than the independent test
when observations are naturally paired, because it eliminates
variability between pairs.
Returns:
Dictionary containing test results
"""
if self.sample2 is None:
raise ValueError("Paired test requires two samples")
if self.n1 != self.n2:
raise ValueError("Paired test requires equal sample sizes")
# Calculate differences for each pair
differences = [self.sample1[i] - self.sample2[i] for i in range(self.n1)]
# Statistics of differences
mean_diff = self._mean(differences)
std_diff = self._std_dev(differences)
standard_error = std_diff / math.sqrt(self.n1)
# T-statistic
if standard_error == 0:
t_statistic = 0
else:
t_statistic = mean_diff / standard_error
# Degrees of freedom
df = self.n1 - 1
# P-value
p_value = self._approximate_t_test_p_value(t_statistic, df)
# Confidence interval
t_critical = self._t_critical_value(df, 0.05)
margin_of_error = t_critical * standard_error
ci_lower = mean_diff - margin_of_error
ci_upper = mean_diff + margin_of_error
return {
'test_type': 'Paired samples t-test',
'n_pairs': self.n1,
'mean_difference': mean_diff,
'std_difference': std_diff,
'std_error': standard_error,
't_statistic': t_statistic,
'degrees_of_freedom': df,
'p_value': p_value,
'significant_at_0.05': p_value < 0.05,
'significant_at_0.01': p_value < 0.01,
'confidence_interval_95': (ci_lower, ci_upper)
}
def _approximate_t_test_p_value(self, t_stat, df):
"""
Approximate two-tailed p-value for t-statistic.
For large df (>30), the t-distribution approximates the
standard normal distribution. For smaller df, we use
conservative approximations based on critical values.
In production code, you would use scipy.stats.t.cdf or
implement the incomplete beta function for exact p-values.
Parameters:
t_stat: The t-statistic
df: Degrees of freedom
Returns:
Approximate p-value
"""
abs_t = abs(t_stat)
if df > 30:
# Use normal approximation for large df
z = abs_t
# Standard normal CDF approximation
p_one_tail = 0.5 * (1 - math.erf(z / math.sqrt(2)))
p_value = 2 * p_one_tail
else:
# Conservative approximation based on critical values
# These are approximate critical values for two-tailed tests
if abs_t > 4:
p_value = 0.001
elif abs_t > 3:
p_value = 0.01
elif abs_t > 2.5:
p_value = 0.02
elif abs_t > 2:
p_value = 0.05
elif abs_t > 1.5:
p_value = 0.15
else:
p_value = 0.30
return min(p_value, 1.0)
def _t_critical_value(self, df, alpha):
"""
Approximate t critical value for two-tailed test.
Returns the t-value such that P(absolute value of T > t) = alpha.
For alpha=0.05 (95% confidence), this is approximately:
- df > 30: t approximately equals 1.96 (approaches normal)
- df = 20: t approximately equals 2.09
- df = 10: t approximately equals 2.23
- df = 5: t approximately equals 2.57
Parameters:
df: Degrees of freedom
alpha: Significance level
Returns:
Approximate critical t-value
"""
if alpha == 0.05:
if df > 30:
return 1.96
elif df > 20:
return 2.09
elif df > 10:
return 2.23
else:
return 2.57
elif alpha == 0.01:
if df > 30:
return 2.58
elif df > 20:
return 2.85
else:
return 3.17
else:
# Default conservative value
return 2.0
def effect_size_cohens_d(self):
"""
Calculate Cohen's d effect size measure.
Cohen's d quantifies the standardized difference between means.
It measures effect size in standard deviation units, making it
independent of the original measurement scale.
Formula: d = (x_bar_1 - x_bar_2) / sp
Where sp is the pooled standard deviation:
sp = square root of [((n_1-1)s_1 squared + (n_2-1)s_2 squared) / (n_1+n_2-2)]
Interpretation (Cohen's guidelines):
- absolute value of d < 0.2: Negligible effect
- 0.2 <= absolute value of d < 0.5: Small effect
- 0.5 <= absolute value of d < 0.8: Medium effect
- absolute value of d >= 0.8: Large effect
Effect size is important because statistical significance
depends on sample size, but effect size measures practical
significance.
Returns:
Float representing Cohen's d
"""
if self.sample2 is None:
raise ValueError("Effect size calculation requires two samples")
mean1 = self._mean(self.sample1)
mean2 = self._mean(self.sample2)
var1 = self._variance(self.sample1)
var2 = self._variance(self.sample2)
# Pooled standard deviation
pooled_var = ((self.n1 - 1) * var1 + (self.n2 - 1) * var2) / (self.n1 + self.n2 - 2)
pooled_std = math.sqrt(pooled_var)
if pooled_std == 0:
return 0
cohens_d = (mean1 - mean2) / pooled_std
return cohens_d
Example usage:
# Example: Comparing response times between two service implementations
# Implementation A response times (milliseconds)
impl_a_times = [
45, 48, 52, 47, 50, 49, 51, 48, 46, 53,
47, 50, 49, 52, 48, 51, 47, 49, 50, 48
]
# Implementation B response times (milliseconds)
impl_b_times = [
38, 42, 40, 39, 43, 41, 38, 44, 40, 42,
39, 41, 40, 43, 38, 42, 41, 39, 40, 42
]
# Perform independent samples t-test
ttest = TTest(impl_a_times, impl_b_times)
results = ttest.independent_samples_test(equal_variance=True)
print("Performance Comparison: Implementation A vs B")
print("=" * 70)
print(f"Test: {results['test_type']}")
print()
print(f"Implementation A:")
print(f" Mean: {results['mean1']:.2f} ms")
print(f" Sample size: {len(impl_a_times)}")
print()
print(f"Implementation B:")
print(f" Mean: {results['mean2']:.2f} ms")
print(f" Sample size: {len(impl_b_times)}")
print()
print(f"Difference (A - B): {results['mean_difference']:.2f} ms")
print(f"95% Confidence Interval: ({results['confidence_interval_95'][0]:.2f}, {results['confidence_interval_95'][1]:.2f}) ms")
print()
print(f"Statistical Test Results:")
print(f" t-statistic: {results['t_statistic']:.4f}")
print(f" Degrees of freedom: {results['degrees_of_freedom']:.1f}")
print(f" p-value: {results['p_value']:.4f}")
print()
print(f"Significance:")
print(f" Significant at alpha=0.05? {results['significant_at_0.05']}")
print(f" Significant at alpha=0.01? {results['significant_at_0.01']}")
print()
# Calculate effect size
effect_size = ttest.effect_size_cohens_d()
print(f"Effect Size:")
print(f" Cohen's d: {effect_size:.4f}")
# Interpret effect size
if abs(effect_size) < 0.2:
interpretation = "negligible"
elif abs(effect_size) < 0.5:
interpretation = "small"
elif abs(effect_size) < 0.8:
interpretation = "medium"
else:
interpretation = "large"
print(f" Interpretation: {interpretation} effect")
print()
# Conclusion
if results['significant_at_0.05']:
print("Conclusion:")
print(f" Implementation B is significantly faster than A (p < 0.05)")
print(f" The difference of {results['mean_difference']:.2f} ms is statistically significant")
print(f" and represents a {interpretation} practical effect.")
else:
print("Conclusion:")
print(f" No significant difference detected (p >= 0.05)")
print(f" Insufficient evidence to conclude implementations differ in performance.")
This comprehensive implementation explains the mathematics and logic of t-tests in detail. The code demonstrates both statistical significance through p-values and practical significance through effect sizes. For software architects, understanding both aspects is crucial for making informed decisions about system changes.
The Chi-Square Test: Analyzing Categorical Data
The chi-square test analyzes categorical data to assess whether observed frequencies differ significantly from expected frequencies. The test is named after the Greek letter chi, which is written as a symbol resembling an X. The test statistic follows a chi-square distribution under the null hypothesis.
The chi-square test for independence examines whether two categorical variables are associated. For example, we might test whether error type is independent of service version, or whether user action depends on system response time category. The test compares observed frequencies in a contingency table against frequencies expected if the variables were truly independent.
The expected frequency for each cell in the contingency table is calculated based on marginal totals. If the row variable and column variable are independent, the expected frequency for a cell equals the row total times the column total divided by the grand total. This represents the frequency we would expect if the row and column proportions were the same across all cells.
The test statistic sums the squared standardized deviations between observed and expected frequencies across all cells. The formula is the sum over all cells of the quantity observed minus expected, squared, divided by expected. This statistic measures how far the observed data deviates from what we would expect under independence.
The chi-square statistic follows a chi-square distribution with degrees of freedom equal to the number of rows minus one, times the number of columns minus one. Large values of the chi-square statistic indicate substantial discrepancy between observation and expectation, providing evidence against independence.
Let‘s implement this:
# Chi-square test implementation with detailed mathematics
import math
class ChiSquareTest:
"""
Implementation of chi-square tests for categorical data.
The chi-square test assesses whether observed frequencies differ
significantly from expected frequencies.
Test statistic:
chi squared = sum of [(O - E) squared / E]
Where:
- O is the observed frequency in each cell
- E is the expected frequency in each cell
- sum means sum over all cells
- chi squared is pronounced "kai squared"
The chi-square statistic follows a chi-square distribution under
the null hypothesis. The chi-square distribution is right-skewed
and takes only positive values.
Types of chi-square tests:
1. Test of independence: Are two categorical variables independent?
2. Goodness of fit: Do observed frequencies match expected proportions?
3. Homogeneity: Do different populations have the same distribution?
"""
def __init__(self, observed):
"""
Initialize with observed frequency table.
Parameters:
observed: 2D list (contingency table) of observed frequencies
Example: [[10, 20], [30, 40]] for 2x2 table
Rows represent one variable, columns represent another
"""
if not observed or not observed[0]:
raise ValueError("Observed frequencies table cannot be empty")
self.observed = observed
self.rows = len(observed)
self.cols = len(observed[0])
# Validate rectangular table
for row in observed:
if len(row) != self.cols:
raise ValueError("All rows must have same number of columns")
# Validate non-negative frequencies
for row in observed:
for freq in row:
if freq < 0:
raise ValueError("Frequencies must be non-negative")
def test_independence(self):
"""
Perform chi-square test of independence.
Tests whether row and column variables are independent.
Null hypothesis: Row and column variables are independent
Alternative hypothesis: Row and column variables are associated
Algorithm:
1. Calculate row totals, column totals, and grand total
2. Calculate expected frequencies:
E[i,j] = (row_total[i] times col_total[j]) / grand_total
3. Calculate chi-square statistic:
chi squared = sum of [(O[i,j] - E[i,j]) squared / E[i,j]]
4. Calculate degrees of freedom: df = (rows - 1) times (cols - 1)
5. Compare chi squared to chi-square distribution with df degrees of freedom
The expected frequency formula comes from independence:
If independent, P(row i and col j) = P(row i) times P(col j)
So E[i,j] = n times P(row i) times P(col j)
= n times (row_total[i]/n) times (col_total[j]/n)
= (row_total[i] times col_total[j]) / n
Returns:
Dictionary containing test results
"""
# Step 1: Calculate marginal totals
row_totals = [sum(row) for row in self.observed]
col_totals = [sum(self.observed[i][j] for i in range(self.rows))
for j in range(self.cols)]
grand_total = sum(row_totals)
if grand_total == 0:
raise ValueError("Total frequency cannot be zero")
# Step 2: Calculate expected frequencies
expected = []
for i in range(self.rows):
expected_row = []
for j in range(self.cols):
# Expected frequency under independence
expected_freq = (row_totals[i] * col_totals[j]) / grand_total
expected_row.append(expected_freq)
expected.append(expected_row)
# Step 3: Calculate chi-square statistic
chi_square = 0.0
for i in range(self.rows):
for j in range(self.cols):
obs = self.observed[i][j]
exp = expected[i][j]
if exp > 0:
# Add (O - E) squared / E for this cell
chi_square += ((obs - exp) ** 2) / exp
elif obs > 0:
# Expected is zero but observed is not: infinite contribution
raise ValueError(f"Expected frequency is zero at cell ({i},{j}) but observed is {obs}")
# Step 4: Degrees of freedom
df = (self.rows - 1) * (self.cols - 1)
# Step 5: Calculate p-value
p_value = self._approximate_chi_square_p_value(chi_square, df)
# Check minimum expected frequency assumption
min_expected = min(min(row) for row in expected)
assumption_met = min_expected >= 5
return {
'test_type': 'Chi-square test of independence',
'chi_square_statistic': chi_square,
'degrees_of_freedom': df,
'p_value': p_value,
'significant_at_0.05': p_value < 0.05,
'significant_at_0.01': p_value < 0.01,
'observed_frequencies': self.observed,
'expected_frequencies': expected,
'row_totals': row_totals,
'col_totals': col_totals,
'grand_total': grand_total,
'min_expected_frequency': min_expected,
'assumption_met': assumption_met,
'assumption_note': 'All expected frequencies should be >= 5'
}
def goodness_of_fit(self, expected_proportions):
"""
Perform chi-square goodness of fit test.
Tests whether observed frequencies match expected proportions.
Null hypothesis: Data follows the specified distribution
Alternative hypothesis: Data does not follow the specified distribution
Algorithm:
1. Convert expected proportions to expected frequencies
2. Calculate chi-square statistic: chi squared = sum of [(O - E) squared / E]
3. Degrees of freedom: df = number of categories - 1
This test is useful for checking whether data follows a
theoretical distribution, such as uniform distribution or
a specified probability model.
Parameters:
expected_proportions: List of expected proportions (must sum to 1)
Returns:
Dictionary containing test results
"""
# Flatten observed frequencies (works for 1D or 2D input)
if self.rows == 1:
observed_flat = self.observed[0]
elif self.cols == 1:
observed_flat = [row[0] for row in self.observed]
else:
raise ValueError("Goodness of fit requires 1D frequency table")
n_categories = len(observed_flat)
if len(expected_proportions) != n_categories:
raise ValueError("Expected proportions must match number of categories")
if abs(sum(expected_proportions) - 1.0) > 0.001:
raise ValueError("Expected proportions must sum to 1")
# Calculate total and expected frequencies
total = sum(observed_flat)
expected_frequencies = [p * total for p in expected_proportions]
# Calculate chi-square statistic
chi_square = 0.0
for obs, exp in zip(observed_flat, expected_frequencies):
if exp > 0:
chi_square += ((obs - exp) ** 2) / exp
elif obs > 0:
raise ValueError("Expected frequency is zero but observed is not")
# Degrees of freedom
df = n_categories - 1
# P-value
p_value = self._approximate_chi_square_p_value(chi_square, df)
# Check assumption
min_expected = min(expected_frequencies)
assumption_met = min_expected >= 5
return {
'test_type': 'Chi-square goodness of fit test',
'chi_square_statistic': chi_square,
'degrees_of_freedom': df,
'p_value': p_value,
'significant_at_0.05': p_value < 0.05,
'significant_at_0.01': p_value < 0.01,
'observed_frequencies': observed_flat,
'expected_frequencies': expected_frequencies,
'min_expected_frequency': min_expected,
'assumption_met': assumption_met
}
def _approximate_chi_square_p_value(self, chi_square, df):
"""
Approximate p-value for chi-square statistic.
The chi-square distribution is the sum of squared standard
normal variables. For df degrees of freedom, it's the sum
of df independent squared standard normals.
Properties:
- Mean = df
- Variance = 2 times df
- Right-skewed, approaches normal as df increases
This function uses conservative approximations based on
critical values from chi-square tables.
In production, use scipy.stats.chi2.sf(chi_square, df) for
exact p-values.
Parameters:
chi_square: The chi-square statistic
df: Degrees of freedom
Returns:
Approximate p-value
"""
# Critical values for common df at alpha = 0.05 and alpha = 0.01
# These are from standard chi-square tables
critical_values_05 = {
1: 3.84, 2: 5.99, 3: 7.81, 4: 9.49, 5: 10.83,
6: 12.59, 7: 14.07, 8: 15.51, 9: 16.92, 10: 18.31
}
critical_values_01 = {
1: 6.63, 2: 9.21, 3: 11.34, 4: 13.28, 5: 15.09,
6: 16.81, 7: 18.48, 8: 20.09, 9: 21.67, 10: 23.21
}
# Get critical values for this df
if df in critical_values_05:
cv_05 = critical_values_05[df]
cv_01 = critical_values_01[df]
else:
# Approximation for large df using normal approximation
# chi squared approximately equals Normal(df, 2 times df) for large df
cv_05 = df + 1.645 * math.sqrt(2 * df)
cv_01 = df + 2.326 * math.sqrt(2 * df)
# Determine approximate p-value
if chi_square < cv_05:
# Not significant at 0.05
if chi_square < df:
return 0.50 # Below mean
else:
return 0.10 # Between mean and 0.05 critical value
elif chi_square < cv_01:
# Significant at 0.05 but not 0.01
return 0.03
else:
# Significant at 0.01
if chi_square > cv_01 * 1.5:
return 0.001
else:
return 0.005
def cramers_v(self):
"""
Calculate Cramér's V effect size measure.
Cramér's V measures the strength of association between
two categorical variables, ranging from 0 (no association)
to 1 (perfect association).
Formula: V = square root of [chi squared / (n times min(rows-1, cols-1))]
Where:
- chi squared is the chi-square statistic
- n is the total sample size
- min(rows-1, cols-1) is the smaller of the two df components
Interpretation:
- V < 0.1: Negligible association
- 0.1 <= V < 0.3: Weak association
- 0.3 <= V < 0.5: Moderate association
- V >= 0.5: Strong association
Returns:
Float representing Cramér's V
"""
results = self.test_independence()
chi_square = results['chi_square_statistic']
n = results['grand_total']
# Minimum of (rows-1) and (cols-1)
min_dim = min(self.rows - 1, self.cols - 1)
if n == 0 or min_dim == 0:
return 0
cramers_v = math.sqrt(chi_square / (n * min_dim))
return cramers_v
Example usage:
# Example: Analyzing error distribution across microservices
# Contingency table: rows = error types, columns = services
# Error types: Timeout, Connection, Validation
# Services: Service A, Service B, Service C
error_distribution = [
[15, 8, 12], # Timeout errors
[22, 18, 25], # Connection errors
[8, 14, 10] # Validation errors
]
chi_test = ChiSquareTest(error_distribution)
results = chi_test.test_independence()
print("Error Distribution Analysis Across Services")
print("=" * 70)
print("Observed Frequencies:")
print(" Service A Service B Service C Total")
print(f"Timeout: {error_distribution[0][0]:3d} {error_distribution[0][1]:3d} {error_distribution[0][2]:3d} {results['row_totals'][0]:3d}")
print(f"Connection: {error_distribution[1][0]:3d} {error_distribution[1][1]:3d} {error_distribution[1][2]:3d} {results['row_totals'][1]:3d}")
print(f"Validation: {error_distribution[2][0]:3d} {error_distribution[2][1]:3d} {error_distribution[2][2]:3d} {results['row_totals'][2]:3d}")
print(f"Total: {results['col_totals'][0]:3d} {results['col_totals'][1]:3d} {results['col_totals'][2]:3d} {results['grand_total']:3d}")
print()
print("Expected Frequencies (if independent):")
print(" Service A Service B Service C")
for i, row in enumerate(results['expected_frequencies']):
error_type = ["Timeout:", "Connection:", "Validation:"][i]
print(f"{error_type:15s} {row[0]:6.2f} {row[1]:6.2f} {row[2]:6.2f}")
print()
print("Statistical Test Results:")
print(f" Chi-square statistic (chi squared): {results['chi_square_statistic']:.4f}")
print(f" Degrees of freedom: {results['degrees_of_freedom']}")
print(f" p-value: {results['p_value']:.4f}")
print()
print("Significance:")
print(f" Significant at alpha=0.05? {results['significant_at_0.05']}")
print(f" Significant at alpha=0.01? {results['significant_at_0.01']}")
print()
# Effect size
effect_size = chi_test.cramers_v()
print(f"Effect Size:")
print(f" Cramér's V: {effect_size:.4f}")
if effect_size < 0.1:
interpretation = "negligible"
elif effect_size < 0.3:
interpretation = "weak"
elif effect_size < 0.5:
interpretation = "moderate"
else:
interpretation = "strong"
print(f" Interpretation: {interpretation} association")
print()
# Assumption check
print("Assumption Check:")
print(f" Minimum expected frequency: {results['min_expected_frequency']:.2f}")
print(f" Assumption met (all >= 5): {results['assumption_met']}")
print()
# Conclusion
if results['significant_at_0.05']:
print("Conclusion:")
print(" Error type and service are NOT independent (p < 0.05)")
print(" Different services exhibit different error patterns")
print(f" The association is {interpretation} in strength")
else:
print("Conclusion:")
print(" No significant association between error type and service (p >= 0.05)")
print(" Error patterns appear consistent across services")
This comprehensive implementation explains the chi-square test mathematics in detail. The code demonstrates both statistical significance and effect size measurement. For software architects, the chi-square test is essential for analyzing categorical patterns in system behavior, such as error distributions, user actions, or deployment outcomes.
PART FOUR: PRACTICAL APPLICATION - A/B TESTING FRAMEWORK
A/B Testing for Architectural Decisions
A/B testing represents one of the most common experimental designs in software architecture, comparing two variants by randomly assigning users or requests to each variant and measuring outcomes. This methodology enables architects to evaluate architectural changes in production environments with real users and workloads.
The A/B testing framework requires defining clear metrics that operationalize success, such as response time, error rate, or resource utilization. These metrics must be measurable, relevant to business or technical objectives, and sensitive to the architectural change being evaluated.
Sample size determination ensures sufficient statistical power to detect meaningful effects. Testing with too few observations may fail to detect genuine improvements, while testing with excessive observations wastes resources. Power analysis balances Type I error, Type II error, effect size, and sample size to determine appropriate experiment duration.
Here is a complete A/B testing framework implementation:
# A/B testing framework for architectural experiments
import math
import random
class ABTest:
"""
Framework for conducting and analyzing A/B tests in software systems.
Supports metric collection, statistical analysis, and decision making.
A/B testing allows software architects to empirically evaluate
architectural changes by comparing two variants (A and B) under
real-world conditions.
Key concepts:
- Variant A (Control): The baseline or current implementation
- Variant B (Treatment): The new or experimental implementation
- Random assignment: Each observation is randomly assigned to A or B
- Statistical analysis: Determine if observed differences are significant
"""
def __init__(self, name, metric_name, metric_unit):
"""
Initialize an A/B test.
Parameters:
name: Descriptive name for the test
metric_name: Name of the metric being measured
metric_unit: Unit of measurement (e.g., 'ms', 'requests/sec')
"""
self.name = name
self.metric_name = metric_name
self.metric_unit = metric_unit
self.variant_a_data = []
self.variant_b_data = []
def record_observation(self, variant, value):
"""
Record an observation for a variant.
Parameters:
variant: 'A' or 'B'
value: Measured metric value
"""
if variant == 'A':
self.variant_a_data.append(value)
elif variant == 'B':
self.variant_b_data.append(value)
else:
raise ValueError("Variant must be 'A' or 'B'")
def assign_variant(self, split_ratio=0.5):
"""
Randomly assign a variant based on split ratio.
This implements random assignment, a cornerstone of
experimental design. Random assignment ensures that
confounding variables are distributed evenly across
variants, allowing causal inference.
Parameters:
split_ratio: Probability of assigning variant A (default 0.5)
Returns:
'A' or 'B'
"""
return 'A' if random.random() < split_ratio else 'B'
def _calculate_stats(self, data):
"""
Calculate basic statistics for a dataset.
Returns:
Dictionary with count, mean, variance, and std dev
"""
n = len(data)
if n == 0:
return None
mean = sum(data) / n
if n > 1:
variance = sum((x - mean) ** 2 for x in data) / (n - 1)
std_dev = math.sqrt(variance)
else:
variance = 0
std_dev = 0
return {
'n': n,
'mean': mean,
'std_dev': std_dev,
'variance': variance
}
def analyze(self, alpha=0.05):
"""
Perform statistical analysis of A/B test results.
This method conducts a two-sample t-test to determine
whether the difference between variants is statistically
significant.
Parameters:
alpha: Significance level for hypothesis testing (default 0.05)
Returns:
Dictionary containing comprehensive analysis results
"""
if len(self.variant_a_data) < 2 or len(self.variant_b_data) < 2:
raise ValueError("Need at least 2 observations per variant for analysis")
stats_a = self._calculate_stats(self.variant_a_data)
stats_b = self._calculate_stats(self.variant_b_data)
# Calculate difference
mean_diff = stats_a['mean'] - stats_b['mean']
percent_change = (mean_diff / stats_a['mean']) * 100 if stats_a['mean'] != 0 else 0
# Pooled standard error for independent samples
se = math.sqrt(stats_a['variance']/stats_a['n'] + stats_b['variance']/stats_b['n'])
# T-statistic
t_stat = mean_diff / se if se > 0 else 0
# Degrees of freedom (Welch-Satterthwaite)
if stats_a['variance'] > 0 and stats_b['variance'] > 0:
numerator = (stats_a['variance']/stats_a['n'] + stats_b['variance']/stats_b['n']) ** 2
denom_a = (stats_a['variance']/stats_a['n'])**2 / (stats_a['n']-1)
denom_b = (stats_b['variance']/stats_b['n'])**2 / (stats_b['n']-1)
df = numerator / (denom_a + denom_b)
else:
df = stats_a['n'] + stats_b['n'] - 2
# Approximate p-value
p_value = self._approximate_t_test_p_value(t_stat, df)
# Confidence interval for difference
# Using t-distribution critical value approximation
t_critical = 2.0 if df > 30 else 2.5 # Simplified
margin_of_error = t_critical * se
ci_lower = mean_diff - margin_of_error
ci_upper = mean_diff + margin_of_error
# Effect size (Cohen's d)
pooled_std = math.sqrt(((stats_a['n']-1)*stats_a['variance'] +
(stats_b['n']-1)*stats_b['variance']) /
(stats_a['n'] + stats_b['n'] - 2))
cohens_d = mean_diff / pooled_std if pooled_std > 0 else 0
return {
'variant_a': stats_a,
'variant_b': stats_b,
'mean_difference': mean_diff,
'percent_change': percent_change,
'standard_error': se,
't_statistic': t_stat,
'degrees_of_freedom': df,
'p_value': p_value,
'significant': p_value < alpha,
'confidence_interval_95': (ci_lower, ci_upper),
'cohens_d': cohens_d,
'alpha': alpha
}
def _approximate_t_test_p_value(self, t_stat, df):
"""
Approximate two-tailed p-value for t-test.
Parameters:
t_stat: The t-statistic
df: Degrees of freedom
Returns:
Approximate p-value
"""
abs_t = abs(t_stat)
if df > 30:
# Normal approximation
z = abs_t
p_one_tail = 0.5 * (1 - math.erf(z / math.sqrt(2)))
return 2 * p_one_tail
else:
# Conservative approximation
if abs_t > 4:
return 0.001
elif abs_t > 3:
return 0.01
elif abs_t > 2:
return 0.05
elif abs_t > 1.5:
return 0.15
else:
return 0.30
def calculate_required_sample_size(self, effect_size, alpha=0.05, power=0.80):
"""
Calculate required sample size per variant for desired power.
Power analysis determines how many observations are needed
to detect an effect of a given size with specified confidence.
This uses a simplified approximation:
n approximately equals (2 times (z_alpha + z_beta) squared) / (effect_size squared)
Where:
- z_alpha is the critical value for significance level alpha
- z_beta is the critical value for power (1 - beta)
- effect_size is Cohen's d
Parameters:
effect_size: Expected Cohen's d effect size
alpha: Significance level (default 0.05)
power: Desired statistical power (default 0.80)
Returns:
Required sample size per variant
"""
# Critical values (approximations)
z_alpha = 1.96 if alpha == 0.05 else 2.58 # Two-tailed
z_beta = 0.84 if power == 0.80 else 1.28 # Power 0.80 or 0.90
# Sample size formula
n = (2 * ((z_alpha + z_beta) ** 2)) / (effect_size ** 2)
return int(math.ceil(n))
def generate_report(self, analysis_results):
"""
Generate human-readable report of A/B test results.
Parameters:
analysis_results: Dictionary from analyze() method
Returns:
Formatted string report
"""
report = []
report.append("=" * 70)
report.append(f"A/B Test Report: {self.name}")
report.append("=" * 70)
report.append("")
report.append(f"Metric: {self.metric_name} ({self.metric_unit})")
report.append("")
report.append("Variant A (Control):")
report.append(f" Sample size: {analysis_results['variant_a']['n']}")
report.append(f" Mean: {analysis_results['variant_a']['mean']:.4f} {self.metric_unit}")
report.append(f" Std Dev: {analysis_results['variant_a']['std_dev']:.4f} {self.metric_unit}")
report.append("")
report.append("Variant B (Treatment):")
report.append(f" Sample size: {analysis_results['variant_b']['n']}")
report.append(f" Mean: {analysis_results['variant_b']['mean']:.4f} {self.metric_unit}")
report.append(f" Std Dev: {analysis_results['variant_b']['std_dev']:.4f} {self.metric_unit}")
report.append("")
report.append("Statistical Analysis:")
report.append(f" Mean difference (A - B): {analysis_results['mean_difference']:.4f} {self.metric_unit}")
report.append(f" Percent change: {analysis_results['percent_change']:.2f}%")
report.append(f" 95% Confidence Interval: ({analysis_results['confidence_interval_95'][0]:.4f}, {analysis_results['confidence_interval_95'][1]:.4f})")
report.append(f" t-statistic: {analysis_results['t_statistic']:.4f}")
report.append(f" p-value: {analysis_results['p_value']:.4f}")
report.append(f" Cohen's d: {analysis_results['cohens_d']:.4f}")
report.append("")
report.append(f"Decision (alpha = {analysis_results['alpha']}):")
if analysis_results['significant']:
report.append(" SIGNIFICANT DIFFERENCE DETECTED")
if analysis_results['mean_difference'] > 0:
report.append(" Variant B performs better (lower is better assumed)")
else:
report.append(" Variant A performs better (lower is better assumed)")
else:
report.append(" NO SIGNIFICANT DIFFERENCE")
report.append(" Insufficient evidence to prefer one variant")
return "\n".join(report)
Example usage:
# Example: A/B testing database connection pooling configuration
# Simulate an A/B test comparing two connection pool sizes
ab_test = ABTest(
name="Database Connection Pool Size Optimization",
metric_name="Query Response Time",
metric_unit="ms"
)
# Simulate data collection
# Variant A: Pool size 10 (control)
# Variant B: Pool size 20 (treatment)
random.seed(42) # For reproducibility
# Simulate 100 requests for each variant
for _ in range(100):
# Variant A: mean 45ms, std dev 8ms
variant_a_time = random.gauss(45, 8)
ab_test.record_observation('A', max(1, variant_a_time))
# Variant B: mean 38ms, std dev 7ms (improved performance)
variant_b_time = random.gauss(38, 7)
ab_test.record_observation('B', max(1, variant_b_time))
# Analyze results
results = ab_test.analyze(alpha=0.05)
# Generate and print report
report = ab_test.generate_report(results)
print(report)
print()
print("=" * 70)
# Power analysis for future tests
print("\nPower Analysis:")
print("To detect a small effect (Cohen's d = 0.3) with 80% power:")
required_n = ab_test.calculate_required_sample_size(effect_size=0.3, alpha=0.05, power=0.80)
print(f" Required sample size per variant: {required_n}")
print("\nTo detect a medium effect (Cohen's d = 0.5) with 80% power:")
required_n = ab_test.calculate_required_sample_size(effect_size=0.5, alpha=0.05, power=0.80)
print(f" Required sample size per variant: {required_n}")
print("\nTo detect a large effect (Cohen's d = 0.8) with 80% power:")
required_n = ab_test.calculate_required_sample_size(effect_size=0.8, alpha=0.05, power=0.80)
print(f" Required sample size per variant: {required_n}")
This A/B testing framework provides software architects with a complete workflow for experimental evaluation of architectural changes. The implementation handles random assignment, data collection, statistical analysis, and interpretation. The framework encapsulates complexity while exposing key parameters like significance level, enabling architects to make informed decisions based on empirical evidence.
CONCLUSION AND SUMMARY
Understanding Statistical Terms and Algorithms
Throughout this tutorial, we have explored fundamental statistical concepts, probability distributions, and hypothesis testing methods. Let me summarize the key terms and algorithms we've covered:
STATISTICAL TERMS EXPLAINED:
- Mean (Average): The sum of all values divided by the count. It represents the balance point of the data.
- Median: The middle value when data is sorted. It divides the dataset into two equal halves.
- Percentile: The value below which a given percentage of observations fall. The 95th percentile is the value below which 95% of data lies.
- Variance: The average squared deviation from the mean. It measures spread but has squared units.
- Standard Deviation: The square root of variance. It measures spread in the original units of the data.
- Coefficient of Variation: The ratio of standard deviation to mean, expressed as a percentage. It enables comparison of variability across different scales.
- Z-Score: The number of standard deviations a value is from the mean. Formula: z = (x - mean) / standard deviation.
- P-Value: The probability of observing data at least as extreme as what we observed, assuming the null hypothesis is true. Small p-values provide evidence against the null hypothesis.
- Alpha (Significance Level): The threshold for rejecting the null hypothesis, typically 0.05 or 0.01. It controls the Type I error rate.
- Degrees of Freedom (df): The number of independent pieces of information available to estimate a parameter. For variance, df = n - 1.
- Confidence Interval: A range of values that likely contains the true population parameter with a specified level of confidence (e.g., 95%).
- Effect Size: A standardized measure of the magnitude of a difference or relationship, independent of sample size. Cohen's d is a common effect size measure.
- Statistical Power: The probability of correctly detecting an effect when it truly exists. Power = 1 - Type II error rate.
- Chi (χ): The Greek letter used to denote the chi-square statistic and distribution.
- Lambda (λ): The rate parameter in exponential and Poisson distributions.
- Mu (μ): The population mean.
- Sigma (σ): The population standard deviation.
KEY ALGORITHMS EXPLAINED:
- Percentile Calculation Algorithm:
- Sort the data
- Calculate position: (percentile / 100) × (n - 1)
- If position is integer, return data[position]
- If fractional, interpolate between adjacent values
- T-Test Algorithm:
- Calculate means and variances for both groups
- Calculate pooled variance (for equal variance assumption)
- Calculate standard error
- Calculate t-statistic: (mean1 - mean2) / standard error
- Determine degrees of freedom
- Compare to t-distribution to get p-value
- Chi-Square Test Algorithm:
- Calculate row and column totals
- Calculate expected frequencies: (row total × column total) / grand total
- Calculate chi-square statistic: Σ[(observed - expected)² / expected]
- Determine degrees of freedom: (rows - 1) × (columns - 1)
- Compare to chi-square distribution to get p-value
- Normal Distribution CDF Algorithm:
- Standardize to z-score: z = (x - μ) / σ
- Apply error function: erf(z / √2)
- Transform: (1 + erf(z / √2)) / 2
- Exponential Distribution Inverse CDF Algorithm:
- Given probability p
- Calculate: x = -ln(1 - p) / λ
- This gives the time value at probability p
PRACTICAL APPLICATIONS FOR SOFTWARE ARCHITECTS:
- Performance Analysis: Use mean, median, and percentiles to characterize response times and identify tail latency issues.
- Capacity Planning: Use probability distributions (Poisson, Normal) to model request arrivals and determine capacity requirements.
- Reliability Engineering: Use exponential distribution to model failure rates and calculate system reliability over time.
- A/B Testing: Use t-tests to rigorously compare architectural alternatives and make data-driven decisions.
- Anomaly Detection: Use statistical thresholds (e.g., 3-sigma rule) to identify unusual system behavior.
- SLA Definition: Use percentiles (P95, P99) to set realistic service level objectives that account for variability.
- Error Analysis: Use chi-square tests to analyze categorical error patterns across services or components.
BEST PRACTICES:
- Always visualize your data before applying statistical tests
- Check assumptions (normality, independence, equal variance)
- Report both statistical significance (p-value) and practical significance (effect size)
- Use confidence intervals to quantify uncertainty
- Consider multiple comparisons when running many tests
- Ensure adequate sample sizes through power analysis
- Document your statistical methodology for reproducibility
This tutorial has provided you with a comprehensive foundation in statistics for software architecture. By understanding these concepts, terms, and algorithms, you can make evidence-based architectural decisions, validate system performance, and communicate results with precision and confidence.
Remember: statistics is a tool for reasoning under uncertainty. It provides powerful methods for extracting insights from data, but human judgment, domain knowledge, and engineering intuition remain essential for interpreting results and making sound architectural decisions.
No comments:
Post a Comment