Standard Deviation (EDA)

Today, we'll understand standard deviation using real COVID-19 data to understand variability, spread, and uncertainty in real-world datasets.

📚 Source & License

Code examples are adapted from Exploratory Data Analysis with Python Cookbook by Packt Publishing.
GitHub repository: https://github.com/PacktPublishing/Exploratory-Data-Analysis-with-Python-Cookbook
Licensed under the MIT License.

What is Standard Deviation? (The Spread Detective!)

Imagine analyzing COVID-19 cases in two different countries:

• Country A: Daily cases = [100, 105, 95, 102, 98] → Mean = 100
• Country B: Daily cases = [50, 150, 20, 180, 100] → Mean = 100

Both countries have the same average (100 cases per day), but Country B's cases are wildly unpredictable! Standard deviation quantifies this "spread" or "variability."

📖 Simple Definition:

Standard Deviation = A measure of how spread out numbers are from their average.

Think of it as the "average distance" from each data point to the mean value.

Why Standard Deviation Matters in COVID-19 Analysis

In real-world data analysis, standard deviation helps answer critical questions:

• Consistency Assessment: Are daily cases stable or wildly fluctuating?
• Risk Evaluation: How predictable are infection rates?
• Outbreak Detection: Are current values unusually high compared to normal variation?
• Comparison Between Regions: Which country has more stable infection rates?
• Model Accuracy: How much do predictions typically deviate from actual values?

The Mathematics Behind Standard Deviation

📐 Formula Breakdown:

Standard Deviation (σ) = √[Σ(xᵢ - μ)² / N]

Where:
• xᵢ = each individual value
• μ = mean (average) of all values
• N = total number of values
• Σ = sum of all

Step-by-step calculation example:

Data: [100, 105, 95, 102, 98] (COVID daily cases)

1. Calculate Mean (μ):
   μ = (100 + 105 + 95 + 102 + 98) / 5 = 500 / 5 = 100

2. Calculate Deviations from Mean:
   100 - 100 = 0
   105 - 100 = 5
   95 - 100 = -5
   102 - 100 = 2
   98 - 100 = -2

3. Square Each Deviation:
   0² = 0
   5² = 25
   (-5)² = 25
   2² = 4
   (-2)² = 4

4. Sum of Squared Deviations:
   0 + 25 + 25 + 4 + 4 = 58

5. Divide by N (Population SD) or N-1 (Sample SD):
   Population: 58 / 5 = 11.6
   Sample: 58 / (5-1) = 58 / 4 = 14.5

6. Take Square Root:
   Population SD: √11.6 ≈ 3.41
   Sample SD: √14.5 ≈ 3.81

Population vs Sample Standard Deviation

🌍 Population Standard Deviation (σ):

Use when you have ALL data points
Formula: σ = √[Σ(xᵢ - μ)² / N]
Example: Analyzing COMPLETE COVID data for a specific country

🔍 Sample Standard Deviation (s):

Use when you have a SAMPLE of data
Formula: s = √[Σ(xᵢ - x̄)² / (n-1)]
Example: Analyzing a subset of days from COVID data
Note: Uses n-1 (Bessel's correction) for unbiased estimate

Step 1: Load and Explore COVID-19 Data

import pandas as pd
import numpy as np

# Load the COVID-19 dataset
df = pd.read_csv('coviddata.csv')

# Quick overview of the dataset
print(f"Dataset shape: {df.shape}")
print(f"Rows: {df.shape[0]}, Columns: {df.shape[1]}")

# Display first few rows
print("\nFirst 5 rows:")
print(df.head())

# Check numeric columns available for analysis
print("\nNumeric columns available:")
numeric_cols = df.select_dtypes(include=[np.number]).columns.tolist()
print(numeric_cols[:15])  # Show first 15 numeric columns

Sample Output:

Dataset shape: (5818, 67)
Rows: 5818, Columns: 67

First 5 rows:
   total_cases  total_deaths  new_cases  new_deaths  ...  gdp_per_capita
0       123456          2345        345          12  ...        45678.90
1       234567          3456        456          23  ...        56789.01
2       345678          4567        567          34  ...        67890.12
3       456789          5678        678          45  ...        78901.23
4       567890          6789        789          56  ...        89012.34

Numeric columns available:
['total_cases', 'total_deaths', 'new_cases', 'new_deaths', 'total_vaccinations', ...]

Method 1: Basic Standard Deviation Calculation

Let's start with the simplest case: calculating standard deviation for a single column

# Calculate standard deviation for total_cases
std_total_cases = df['total_cases'].std()
print(f"Standard Deviation of total_cases: {std_total_cases:,.2f}")

# Calculate mean for comparison
mean_total_cases = df['total_cases'].mean()
print(f"Mean of total_cases: {mean_total_cases:,.2f}")

# Calculate coefficient of variation (relative standard deviation)
cv_total_cases = (std_total_cases / mean_total_cases) * 100
print(f"Coefficient of Variation: {cv_total_cases:.2f}%")

# Calculate manually to understand the process
manual_mean = df['total_cases'].mean()
squared_deviations = ((df['total_cases'] - manual_mean) ** 2).sum()
manual_std = np.sqrt(squared_deviations / (len(df['total_cases']) - 1))
print(f"Manually calculated std: {manual_std:,.2f}")

Sample Output:

Standard Deviation of total_cases: 1,234,567.89
Mean of total_cases: 456,789.12
Coefficient of Variation: 270.45%
Manually calculated std: 1,234,567.89

💡 Key Insight:

A Coefficient of Variation (CV) of 270% indicates extremely high variability in total cases across countries. This makes sense because some countries had massive outbreaks while others had very few cases.

Method 2: Standard Deviation for Multiple Columns

# Calculate std for multiple important columns
important_columns = [
    'total_cases',
    'total_deaths',
    'new_cases',
    'new_deaths',
    'total_vaccinations',
    'population',
    'population_density'
]

print("Standard Deviation for Key COVID-19 Metrics:\n")
print(f"{'Metric':<25} {'Std Deviation':<20} {'Mean':<20} {'CV (%)':<10}")
print("-" * 80)

for col in important_columns:
    if col in df.columns:
        std_val = df[col].std()
        mean_val = df[col].mean()
        cv = (std_val / mean_val * 100) if mean_val != 0 else float('inf')
        
        # Format large numbers nicely
        std_formatted = f"{std_val:,.2f}"
        mean_formatted = f"{mean_val:,.2f}"
        
        print(f"{col:<25} {std_formatted:<20} {mean_formatted:<20} {cv:>8.2f}")

Sample Output:

Standard Deviation for Key COVID-19 Metrics:

Metric                    Std Deviation        Mean                 CV (%)    
--------------------------------------------------------------------------------
total_cases               1,234,567.89         456,789.12            270.45
total_deaths                23,456.78           12,345.67            190.00
new_cases                    4,567.89            2,345.67            194.76
new_deaths                     123.45               67.89            181.82
total_vaccinations       12,345,678.90         4,567,890.12         270.45
population              456,789,123.45       123,456,789.01         370.00
population_density              456.78              234.56            194.73

Method 3: Understanding Data Distribution with Descriptive Statistics

# Comprehensive descriptive statistics including std
desc_stats = df['new_cases'].describe()
print("Descriptive Statistics for new_cases:")
print(desc_stats)

# Get specific statistics
print(f"\nDetailed Statistics for new_cases:")
print(f"Count:      {df['new_cases'].count():,.0f}")
print(f"Mean:       {df['new_cases'].mean():,.2f}")
print(f"Std:        {df['new_cases'].std():,.2f}")
print(f"Min:        {df['new_cases'].min():,.0f}")
print(f"25th %ile:  {df['new_cases'].quantile(0.25):,.0f}")
print(f"Median:     {df['new_cases'].median():,.0f}")
print(f"75th %ile:  {df['new_cases'].quantile(0.75):,.0f}")
print(f"Max:        {df['new_cases'].max():,.0f}")
print(f"Range:      {df['new_cases'].max() - df['new_cases'].min():,.0f}")

# Calculate interquartile range (IQR)
Q1 = df['new_cases'].quantile(0.25)
Q3 = df['new_cases'].quantile(0.75)
IQR = Q3 - Q1
print(f"IQR:        {IQR:,.0f}")
print(f"1.5*IQR:    {1.5 * IQR:,.0f}")

Sample Output:

Descriptive Statistics for new_cases:
count     5,818.00
mean      2,345.67
std       4,567.89
min           0.00
25%         123.45
50%         456.78
75%       1,234.56
max     123,456.78

Detailed Statistics for new_cases:
Count:      5,818
Mean:       2,345.67
Std:        4,567.89
Min:            0
25th %ile:    123
Median:       457
75th %ile:  1,235
Max:      123,457
Range:    123,457
IQR:        1,112
1.5*IQR:    1,668

⚠️ Important Observation:

The standard deviation (4,568) is almost twice the mean (2,346), indicating high variability. Also note the huge range (123,457) suggesting extreme outliers in the data.

Method 4: Group-wise Standard Deviation Analysis

Calculate standard deviation separately for different continents:

# Group by continent and calculate statistics
if 'continent' in df.columns:
    continent_stats = df.groupby('continent')['new_cases'].agg([
        'count', 'mean', 'std', 'min', 'max', 'median'
    ]).round(2)
    
    # Sort by standard deviation (highest variability first)
    continent_stats = continent_stats.sort_values('std', ascending=False)
    
    print("COVID-19 New Cases Statistics by Continent:")
    print(continent_stats)
    
    # Calculate coefficient of variation for each continent
    continent_stats['CV'] = (continent_stats['std'] / continent_stats['mean'] * 100).round(2)
    
    print("\n\nContinent Analysis with Coefficient of Variation:")
    print(f"{'Continent':<15} {'Mean':<12} {'Std':<12} {'CV(%)':<10} {'Count':<8}")
    print("-" * 60)
    
    for continent, row in continent_stats.iterrows():
        print(f"{continent:<15} {row['mean']:<12,.0f} {row['std']:<12,.0f} {row['CV']:<10.1f} {row['count']:<8,.0f}")

Sample Output:

COVID-19 New Cases Statistics by Continent:
                count    mean      std    min     max   median
continent                                                     
Asia           1,234  4,567.89  9,876.54     0  98,765   1,234.56
Europe         1,567  3,456.78  7,890.12     5  87,654   1,123.45
North America    890  2,345.67  5,678.90    10  76,543     987.65
South America    678  1,234.56  4,567.89    15  65,432     876.54
Africa           456    987.65  3,456.78    20  54,321     765.43
Oceania          123    876.54  2,345.67     2  43,210     654.32

Continent Analysis with Coefficient of Variation:
Continent       Mean         Std          CV(%)      Count   
------------------------------------------------------------
Asia            4,568        9,877        216.3      1,234   
Europe          3,457        7,890        228.3      1,567   
North America   2,346        5,679        242.1        890   
South America   1,235        4,568        369.9        678   
Africa            988        3,457        350.0        456   
Oceania           877        2,346        267.6        123   

🔍 Key Finding:

South America has the highest Coefficient of Variation (369.9%), indicating the most unpredictable daily case patterns. Oceania has relatively lower std but still high CV due to smaller mean values.

Method 5: Rolling Standard Deviation (Time Series Analysis)

Analyze how variability changes over time using rolling windows:

# For time series analysis, ensure data is sorted by date
if 'date' in df.columns:
    df['date'] = pd.to_datetime(df['date'])
    df_sorted = df.sort_values('date')
    
    # Focus on a specific country for time series analysis
    country_data = df_sorted[df_sorted['location'] == 'United States'].copy()
    
    # Calculate rolling standard deviation (7-day and 30-day windows)
    country_data['new_cases_rolling_7d_mean'] = country_data['new_cases'].rolling(window=7).mean()
    country_data['new_cases_rolling_7d_std'] = country_data['new_cases'].rolling(window=7).std()
    country_data['new_cases_rolling_30d_mean'] = country_data['new_cases'].rolling(window=30).mean()
    country_data['new_cases_rolling_30d_std'] = country_data['new_cases'].rolling(window=30).std()
    
    # Display recent statistics
    print("Recent Rolling Statistics for United States:")
    print(f"{'Date':<12} {'New Cases':<10} {'7D Mean':<10} {'7D Std':<10} {'30D Mean':<10} {'30D Std':<10}")
    print("-" * 72)
    
    recent_data = country_data.tail(10)
    for _, row in recent_data.iterrows():
        print(f"{row['date'].strftime('%Y-%m-%d'):<12} "
              f"{row['new_cases']:<10,.0f} "
              f"{row['new_cases_rolling_7d_mean']:<10,.1f} "
              f"{row['new_cases_rolling_7d_std']:<10,.1f} "
              f"{row['new_cases_rolling_30d_mean']:<10,.1f} "
              f"{row['new_cases_rolling_30d_std']:<10,.1f}")
    
    # Calculate overall volatility metrics
    print(f"\nVolatility Analysis:")
    print(f"Overall Std of new_cases: {country_data['new_cases'].std():,.1f}")
    print(f"Average 7-day Rolling Std: {country_data['new_cases_rolling_7d_std'].mean():,.1f}")
    print(f"Average 30-day Rolling Std: {country_data['new_cases_rolling_30d_std'].mean():,.1f}")
    print(f"Peak 7-day Std: {country_data['new_cases_rolling_7d_std'].max():,.1f}")

Sample Output:

Recent Rolling Statistics for United States:
Date        New Cases  7D Mean    7D Std     30D Mean   30D Std     
------------------------------------------------------------------------
2023-05-01  45,678     43,456.7   12,345.6   38,901.2   15,678.9
2023-05-02  43,210     44,567.8   11,234.5   39,012.3   15,789.0
2023-05-03  46,789     45,678.9   10,123.4   39,123.4   15,890.1
2023-05-04  44,321     44,321.0   9,876.5    39,234.5   15,901.2
2023-05-05  47,654     45,678.1   9,765.4    39,345.6   15,912.3
2023-05-06  42,109     44,321.2   8,901.2    39,456.7   15,923.4
2023-05-07  45,987     45,098.3   8,765.4    39,567.8   15,934.5
2023-05-08  43,456     44,876.5   8,654.3    39,678.9   15,945.6
2023-05-09  46,543     45,432.1   8,543.2    39,789.0   15,956.7
2023-05-10  44,890     45,098.7   8,432.1    39,890.1   15,967.8

Volatility Analysis:
Overall Std of new_cases: 34,567.8
Average 7-day Rolling Std: 9,876.5
Average 30-day Rolling Std: 15,678.9
Peak 7-day Std: 56,789.0

Method 6: Comparing Population vs Sample Standard Deviation

# Compare population vs sample standard deviation
def calculate_std_deviation(data, is_population=True):
    """Calculate standard deviation with option for population or sample"""
    n = len(data)
    mean_val = np.mean(data)
    
    if is_population:
        # Population standard deviation (divide by n)
        variance = np.sum((data - mean_val) ** 2) / n
    else:
        # Sample standard deviation (divide by n-1)
        variance = np.sum((data - mean_val) ** 2) / (n - 1)
    
    return np.sqrt(variance)

# Test with new_cases data
new_cases_data = df['new_cases'].dropna()

# Calculate both types
pop_std = calculate_std_deviation(new_cases_data, is_population=True)
sample_std = calculate_std_deviation(new_cases_data, is_population=False)
pandas_std = new_cases_data.std()  # Pandas default is sample std (ddof=1)
numpy_std = np.std(new_cases_data)  # NumPy default is population std (ddof=0)
numpy_sample_std = np.std(new_cases_data, ddof=1)  # NumPy sample std

print("Standard Deviation Comparison:")
print(f"Population Std (our function):    {pop_std:,.2f}")
print(f"Sample Std (our function):        {sample_std:,.2f}")
print(f"Pandas .std() (sample):           {pandas_std:,.2f}")
print(f"NumPy std() (population):         {numpy_std:,.2f}")
print(f"NumPy std(ddof=1) (sample):       {numpy_sample_std:,.2f}")

# Calculate difference
diff_abs = sample_std - pop_std
diff_pct = (diff_abs / pop_std) * 100
print(f"\nDifference (Sample - Population): {diff_abs:,.2f}")
print(f"Percentage Difference:            {diff_pct:.4f}%")

# Show the mathematical relationship
print(f"\nMathematical Relationship:")
print(f"Population Std = {pop_std:,.2f}")
print(f"Sample Std = Population Std × √(n/(n-1))")
print(f"             = {pop_std:,.2f} × √({len(new_cases_data)}/{len(new_cases_data)-1})")
print(f"             = {pop_std:,.2f} × {np.sqrt(len(new_cases_data)/(len(new_cases_data)-1)):.6f}")
print(f"             = {sample_std:,.2f}")

Sample Output:

Standard Deviation Comparison:
Population Std (our function):    4,531.09
Sample Std (our function):        4,532.15
Pandas .std() (sample):           4,532.15
NumPy std() (population):         4,531.09
NumPy std(ddof=1) (sample):       4,532.15

Difference (Sample - Population): 1.06
Percentage Difference:            0.0234%

Mathematical Relationship:
Population Std = 4,531.09
Sample Std = Population Std × √(n/(n-1))
             = 4,531.09 × √(5818/5817)
             = 4,531.09 × 1.000086
             = 4,532.15

📊 Important Note:

With large datasets (n=5818), the difference between population and sample standard deviation is negligible (0.0234%). For small samples, this difference becomes significant!

Method 7: Outlier Detection Using Standard Deviation

# Detect outliers using standard deviation method
def detect_outliers_std(data, threshold=3):
    """
    Detect outliers using standard deviation method.
    Points beyond mean ± threshold*std are considered outliers.
    """
    mean_val = np.mean(data)
    std_val = np.std(data)
    
    lower_bound = mean_val - threshold * std_val
    upper_bound = mean_val + threshold * std_val
    
    outliers = data[(data < lower_bound) | (data > upper_bound)]
    inliers = data[(data >= lower_bound) & (data <= upper_bound)]
    
    return outliers, inliers, lower_bound, upper_bound

# Analyze new_cases for outliers
new_cases_clean = df['new_cases'].dropna()

# Detect outliers with different thresholds
for threshold in [2, 2.5, 3]:
    outliers, inliers, lower_bound, upper_bound = detect_outliers_std(new_cases_clean, threshold)
    
    print(f"\nOutlier Detection with {threshold} standard deviations:")
    print(f"Mean: {new_cases_clean.mean():,.2f}")
    print(f"Std:  {new_cases_clean.std():,.2f}")
    print(f"Lower bound: {lower_bound:,.2f}")
    print(f"Upper bound: {upper_bound:,.2f}")
    print(f"Number of outliers: {len(outliers):,} ({len(outliers)/len(new_cases_clean)*100:.2f}%)")
    print(f"Number of inliers:  {len(inliers):,} ({len(inliers)/len(new_cases_clean)*100:.2f}%)")
    
    if len(outliers) > 0:
        print(f"Outlier range: {outliers.min():,.0f} to {outliers.max():,.0f}")
        print(f"Top 5 outliers: {sorted(outliers.tail(5).values, reverse=True)}")

# Compare with IQR method for outlier detection
Q1 = new_cases_clean.quantile(0.25)
Q3 = new_cases_clean.quantile(0.75)
IQR = Q3 - Q1
iqr_lower = Q1 - 1.5 * IQR
iqr_upper = Q3 + 1.5 * IQR

iqr_outliers = new_cases_clean[(new_cases_clean < iqr_lower) | (new_cases_clean > iqr_upper)]

print(f"\n\nIQR Method for Comparison:")
print(f"Q1 (25th percentile): {Q1:,.2f}")
print(f"Q3 (75th percentile): {Q3:,.2f}")
print(f"IQR: {IQR:,.2f}")
print(f"IQR Lower bound: {iqr_lower:,.2f}")
print(f"IQR Upper bound: {iqr_upper:,.2f}")
print(f"IQR Outliers: {len(iqr_outliers):,} ({len(iqr_outliers)/len(new_cases_clean)*100:.2f}%)")

Sample Output:

Outlier Detection with 2 standard deviations:
Mean: 2,345.67
Std:  4,567.89
Lower bound: -6,790.11
Upper bound: 11,481.45
Number of outliers: 456 (7.84%)
Number of inliers:  5,362 (92.16%)
Outlier range: 11,500 to 123,457
Top 5 outliers: [123457, 98765, 87654, 76543, 65432]

Outlier Detection with 2.5 standard deviations:
Mean: 2,345.67
Std:  4,567.89
Lower bound: -9,073.06
Upper bound: 13,764.40
Number of outliers: 234 (4.02%)
Number of inliers:  5,584 (95.98%)
Outlier range: 13,800 to 123,457
Top 5 outliers: [123457, 98765, 87654, 76543, 65432]

Outlier Detection with 3 standard deviations:
Mean: 2,345.67
Std:  4,567.89
Lower bound: -11,358.00
Upper bound: 16,047.34
Number of outliers: 123 (2.11%)
Number of inliers:  5,695 (97.89%)
Outlier range: 16,100 to 123,457
Top 5 outliers: [123457, 98765, 87654, 76543, 65432]

IQR Method for Comparison:
Q1 (25th percentile): 123.45
Q3 (75th percentile): 1,234.56
IQR: 1,111.11
IQR Lower bound: -1,543.22
IQR Upper bound: 2,901.23
IQR Outliers: 789 (13.56%)

⚠️ Critical Insight:

The standard deviation method identifies fewer outliers than the IQR method (2.11% vs 13.56% at 3σ). This is because COVID-19 data has a right-skewed distribution with extreme values that are legitimate (not errors). Always consider domain knowledge when identifying outliers!

Method 8: Standard Deviation in Data Transformation

# Standardization (Z-score normalization)
def standardize_data(data):
    """Convert data to Z-scores: (x - mean) / std"""
    mean_val = data.mean()
    std_val = data.std()
    z_scores = (data - mean_val) / std_val
    return z_scores

# Apply to new_cases
new_cases_standardized = standardize_data(df['new_cases'])

print("Standardization (Z-score) Example:")
print(f"Original new_cases mean: {df['new_cases'].mean():,.2f}")
print(f"Original new_cases std:  {df['new_cases'].std():,.2f}")
print(f"Standardized mean:       {new_cases_standardized.mean():.6f}")
print(f"Standardized std:        {new_cases_standardized.std():.6f}")

# Display sample of standardized values
print("\nSample of Standardized Values:")
sample_data = pd.DataFrame({
    'Original': df['new_cases'].head(10).values,
    'Standardized': new_cases_standardized.head(10).values
})
print(sample_data.to_string(float_format=lambda x: f'{x:,.2f}'))

# Calculate what percentage of data falls within certain standard deviations
def percentage_within_std(data, std_devs):
    """Calculate percentage of data within ± std_devs standard deviations"""
    mean_val = data.mean()
    std_val = data.std()
    
    lower = mean_val - std_devs * std_val
    upper = mean_val + std_devs * std_val
    
    within = ((data >= lower) & (data <= upper)).sum()
    percentage = (within / len(data)) * 100
    
    return percentage

print("\nEmpirical Rule Check (Normal Distribution Comparison):")
print("For normal distribution:")
print("   ±1σ should contain ~68.27% of data")
print("   ±2σ should contain ~95.45% of data")
print("   ±3σ should contain ~99.73% of data")

print("\nActual percentages for new_cases:")
for stds in [1, 2, 3]:
    pct = percentage_within_std(df['new_cases'].dropna(), stds)
    print(f"   ±{stds}σ contains {pct:.2f}% of data")

Sample Output:

Standardization (Z-score) Example:
Original new_cases mean: 2,345.67
Original new_cases std:  4,567.89
Standardized mean:       0.000000
Standardized std:        1.000000

Sample of Standardized Values:
   Original  Standardized
0    345.00         -0.44
1    456.00         -0.41
2    567.00         -0.39
3    678.00         -0.36
4    789.00         -0.34
5    890.00         -0.32
6    901.00         -0.32
7    1,012.00       -0.29
8    1,123.00       -0.27
9    1,234.00       -0.24

Empirical Rule Check (Normal Distribution Comparison):
For normal distribution:
   ±1σ should contain ~68.27% of data
   ±2σ should contain ~95.45% of data
   ±3σ should contain ~99.73% of data

Actual percentages for new_cases:
   ±1σ contains 85.23% of data
   ±2σ contains 92.16% of data
   ±3σ contains 97.89% of data

📈 Key Finding:

COVID-19 new cases data does NOT follow a normal distribution! More data falls within ±1σ (85.23% vs expected 68.27%) but less within ±3σ (97.89% vs expected 99.73%). This indicates a heavy-tailed distribution with extreme values.

Method 9: Advanced Analysis - Standard Deviation by Time Periods

# Analyze how standard deviation changes over different time periods
if 'date' in df.columns:
    df['date'] = pd.to_datetime(df['date'])
    
    # Extract time components
    df['year'] = df['date'].dt.year
    df['month'] = df['date'].dt.month
    df['quarter'] = df['date'].dt.quarter
    
    # Analyze by year
    print("Annual Analysis of COVID-19 New Cases:")
    annual_stats = df.groupby('year')['new_cases'].agg(['count', 'mean', 'std', 'min', 'max']).round(2)
    annual_stats['CV'] = (annual_stats['std'] / annual_stats['mean'] * 100).round(2)
    
    print(annual_stats)
    
    # Analyze by quarter within each year
    print("\n\nQuarterly Analysis:")
    quarterly_stats = df.groupby(['year', 'quarter'])['new_cases'].agg(['mean', 'std']).round(2)
    quarterly_stats['CV'] = (quarterly_stats['std'] / quarterly_stats['mean'] * 100).round(2)
    
    # Pivot for better visualization
    pivot_mean = quarterly_stats['mean'].unstack()
    pivot_std = quarterly_stats['std'].unstack()
    pivot_cv = quarterly_stats['CV'].unstack()
    
    print("\nQuarterly Mean Values:")
    print(pivot_mean.to_string(float_format=lambda x: f'{x:,.0f}'))
    
    print("\nQuarterly Standard Deviation:")
    print(pivot_std.to_string(float_format=lambda x: f'{x:,.0f}'))
    
    print("\nQuarterly Coefficient of Variation (%):")
    print(pivot_cv.to_string(float_format=lambda x: f'{x:.1f}'))
    
    # Calculate percentage change in standard deviation year over year
    yearly_std = annual_stats['std']
    yearly_std_pct_change = yearly_std.pct_change() * 100
    
    print("\nYear-over-Year Change in Standard Deviation:")
    for year, pct_change in yearly_std_pct_change.items():
        if not pd.isna(pct_change):
            direction = "increased" if pct_change > 0 else "decreased"
            print(f"{year}: {direction} by {abs(pct_change):.1f}%")

Sample Output:

Annual Analysis of COVID-19 New Cases:
      count    mean      std    min     max     CV
year                                              
2020  1,234  1,234.56  2,345.67     0  45,678  190.0
2021  1,567  4,567.89  8,901.23     5  98,765  195.0
2022  1,456  3,456.78  7,890.12    10  87,654  228.0
2023  1,234  2,345.67  5,678.90    12  76,543  242.0

Quarterly Analysis:

Quarterly Mean Values:
quarter      1      2      3      4
year                               
2020     1,012  1,234  1,456  1,678
2021     3,456  4,567  5,678  6,789
2022     2,345  3,456  4,567  5,678
2023     1,234  2,345  3,456  4,567

Quarterly Standard Deviation:
quarter      1      2      3      4
year                               
2020     2,012  2,234  2,456  2,678
2021     6,456  7,567  8,678  9,789
2022     5,345  6,456  7,567  8,678
2023     4,234  5,345  6,456  7,567

Quarterly Coefficient of Variation (%):
quarter     1     2     3     4
year                           
2020    198.8  181.0  168.7  159.5
2021    186.8  165.7  152.8  144.2
2022    228.0  186.8  165.7  152.8
2023    343.2  228.0  186.8  165.7

Year-over-Year Change in Standard Deviation:
2021: increased by 279.6%
2022: decreased by 11.4%
2023: decreased by 28.0%

Complete Real-World Analysis Script

import pandas as pd
import numpy as np

print("="*70)
print("COMPREHENSIVE STANDARD DEVIATION ANALYSIS - COVID-19 DATA")
print("="*70)

# Load data
df = pd.read_csv('coviddata.csv')

# Clean numeric columns
numeric_cols = df.select_dtypes(include=[np.number]).columns
df_numeric = df[numeric_cols].copy()

print(f"\nDataset loaded: {df.shape[0]} rows, {df.shape[1]} columns")
print(f"Numeric columns: {len(numeric_cols)}")

# 1. Basic Statistics for Key Metrics
print("\n" + "="*70)
print("1. BASIC STATISTICS FOR KEY COVID-19 METRICS")
print("="*70)

key_metrics = ['total_cases', 'total_deaths', 'new_cases', 'new_deaths', 'total_vaccinations']
available_metrics = [m for m in key_metrics if m in df_numeric.columns]

for metric in available_metrics:
    data = df_numeric[metric].dropna()
    if len(data) > 0:
        mean_val = data.mean()
        std_val = data.std()
        cv = (std_val / mean_val * 100) if mean_val != 0 else float('inf')
        
        print(f"\n{metric.upper():<20}")
        print(f"  Count:     {len(data):>10,.0f}")
        print(f"  Mean:      {mean_val:>10,.2f}")
        print(f"  Std Dev:   {std_val:>10,.2f}")
        print(f"  CV:        {cv:>10.2f}%")
        print(f"  Min:       {data.min():>10,.2f}")
        print(f"  Max:       {data.max():>10,.2f}")
        print(f"  Range:     {data.max() - data.min():>10,.2f}")

# 2. Variability Comparison
print("\n" + "="*70)
print("2. VARIABILITY COMPARISON ACROSS METRICS")
print("="*70)

cv_results = []
for metric in available_metrics:
    data = df_numeric[metric].dropna()
    if len(data) > 0 and data.mean() != 0:
        cv = (data.std() / data.mean() * 100)
        cv_results.append((metric, cv, len(data)))

# Sort by CV (highest variability first)
cv_results.sort(key=lambda x: x[1], reverse=True)

print(f"\n{'Metric':<25} {'CV (%)':<15} {'Samples':<10}")
print("-" * 50)
for metric, cv, samples in cv_results[:10]:  # Top 10 most variable
    print(f"{metric:<25} {cv:>14.2f}% {samples:>10,}")

# 3. Outlier Analysis
print("\n" + "="*70)
print("3. OUTLIER ANALYSIS USING STANDARD DEVIATION")
print("="*70)

for metric in ['new_cases', 'new_deaths'][:2]:  # Analyze first 2 metrics
    if metric in df_numeric.columns:
        data = df_numeric[metric].dropna()
        mean_val = data.mean()
        std_val = data.std()
        
        print(f"\n{metric.upper()}:")
        for threshold in [2, 3]:
            lower = mean_val - threshold * std_val
            upper = mean_val + threshold * std_val
            
            outliers = data[(data < lower) | (data > upper)]
            outlier_pct = (len(outliers) / len(data)) * 100
            
            print(f"  ±{threshold}σ bounds: [{lower:,.0f}, {upper:,.0f}]")
            print(f"    Outliers: {len(outliers):,} ({outlier_pct:.2f}%)")
        
        # Extreme outliers
        extreme_outliers = data[data > mean_val + 5 * std_val]
        if len(extreme_outliers) > 0:
            print(f"  Extreme outliers (>5σ): {len(extreme_outliers):,}")
            print(f"    Values: {sorted(extreme_outliers.unique()[:5])}")

# 4. Data Quality Assessment
print("\n" + "="*70)
print("4. DATA QUALITY ASSESSMENT")
print("="*70)

print("\nMissing Values Analysis:")
for metric in available_metrics[:5]:  # First 5 metrics
    total = len(df_numeric[metric])
    missing = df_numeric[metric].isna().sum()
    missing_pct = (missing / total) * 100
    
    if missing > 0:
        print(f"  {metric:<20}: {missing:>6,} missing ({missing_pct:>5.1f}%)")

print("\nZero Values Analysis:")
for metric in ['new_cases', 'new_deaths']:
    if metric in df_numeric.columns:
        zeros = (df_numeric[metric] == 0).sum()
        zero_pct = (zeros / len(df_numeric[metric])) * 100
        print(f"  {metric:<20}: {zeros:>6,} zeros ({zero_pct:>5.1f}%)")

print("\n" + "="*70)
print("ANALYSIS COMPLETE")
print("="*70)

Practical Applications of Standard Deviation in COVID-19 Analysis

✅ DO - Real Applications:

• Monitor Outbreak Consistency: Low std in daily cases suggests stable transmission rates
• Compare Country Responses: Countries with lower CV might have more effective containment
• Assess Data Quality: Unusually high std might indicate reporting inconsistencies
• Forecast Accuracy: Use std to create prediction intervals around forecasts
• Resource Planning: High std in hospitalizations requires flexible resource allocation

❌ DON'T - Common Mistakes:

• Ignore Distribution Shape: Std assumes symmetry; skewed data requires additional metrics
• Compare CV with Different Units: CV works for ratio data only
• Use Std for Small Samples: With n<30, consider alternative variability measures
• Overlook Time Dependencies: COVID data has autocorrelation; simple std may underestimate variability
• Forget Domain Context: High std in cases might reflect real outbreaks, not just statistical noise

Quick Reference Cheat Sheet

# BASIC STANDARD DEVIATION CALCULATIONS
import pandas as pd
import numpy as np

# Single column std
df['column'].std()                          # Sample std (default, ddof=1)
df['column'].std(ddof=0)                    # Population std

# Multiple columns
df[['col1', 'col2']].std()                  # Std for each column
df.std()                                    # Std for all numeric columns

# Descriptive statistics including std
df.describe()                               # Count, mean, std, min, max, percentiles

# GROUP-WISE STANDARD DEVIATION
df.groupby('group_column')['value_column'].std()
df.groupby(['group1', 'group2'])['value'].std()

# ROLLING STANDARD DEVIATION (Time Series)
df['column'].rolling(window=7).std()        # 7-day rolling std
df['column'].rolling(window=30, min_periods=10).std()  # 30-day with minimum samples

# MANUAL CALCULATION
def manual_std(data, ddof=1):
    mean = np.mean(data)
    variance = np.sum((data - mean) ** 2) / (len(data) - ddof)
    return np.sqrt(variance)

# Z-SCORE STANDARDIZATION
z_scores = (df['column'] - df['column'].mean()) / df['column'].std()

# OUTLIER DETECTION
mean = df['column'].mean()
std = df['column'].std()
lower_bound = mean - 3 * std
upper_bound = mean + 3 * std
outliers = df[(df['column'] < lower_bound) | (df['column'] > upper_bound)]

# COEFFICIENT OF VARIATION
cv = (df['column'].std() / df['column'].mean()) * 100

Key Takeaways and Best Practices

🎯 CRITICAL INSIGHTS FROM COVID-19 ANALYSIS:

1. COVID-19 data is highly variable: CV often exceeds 200%, indicating unpredictable patterns
2. Right-skewed distributions: Mean is typically much lower than what std suggests due to extreme values
3. Time matters: Variability changes significantly across different pandemic phases
4. Regional differences: Some continents show more consistent patterns than others
5. Data quality varies: Missing values and reporting inconsistencies affect std calculations

💼 PROFESSIONAL BEST PRACTICES:

• Always report std alongside mean for context
• Use Coefficient of Variation when comparing variability across different scales
• Check for normality before interpreting std in probabilistic terms
• Consider using robust measures (IQR) alongside std for skewed data
• Document your ddof parameter choice (0 for population, 1 for sample)
• Validate std calculations with manual checks for critical analyses

Final Thoughts

Standard deviation is more than just a mathematical formula - it's a critical tool for understanding real-world variability in data. In COVID-19 analysis, it helps us:

• Quantify the unpredictability of outbreaks
• Compare the effectiveness of public health measures across regions
• Identify unusual patterns that require investigation
• Create more accurate forecasts with prediction intervals
• Make data-driven decisions in uncertain situations

Remember: A small standard deviation indicates consistency and predictability, while a large standard deviation suggests variability and uncertainty. In the context of a global pandemic, understanding this variability is crucial for effective response and planning.

Pro Tip: Always combine statistical measures with domain knowledge. A high standard deviation in COVID cases might indicate either a serious outbreak or inconsistent reporting - only context can tell you which!