How to Perform Basic Statistical Analysis Using NumPy in Python

pip install numpy

!pip install numpy

import numpy as np

# Creating a NumPy array from a list
my_data = np.array([1, 2, 3, 4, 5])
print(my_data)
print(type(my_data))

[1 2 3 4 5]
<class 'numpy.ndarray'>

# Create an array of zeros
zeros_array = np.zeros(5)  # Creates [0., 0., 0., 0., 0.]
print("Zeros array:", zeros_array)

# Create an array of ones
ones_array = np.ones(5)    # Creates [1., 1., 1., 1., 1.]
print("Ones array:", ones_array)

# Create a range of numbers
range_array = np.arange(0, 10, 2)  # Start at 0, go to 10 (exclusive), step by 2
print("Range array:", range_array)  # [0, 2, 4, 6, 8]

# Create evenly spaced numbers
linspace_array = np.linspace(0, 1, 5)  # 5 numbers evenly spaced between 0 and 1
print("Linspace array:", linspace_array)  # [0., 0.25, 0.5, 0.75, 1.]

# Create a simple 1D array
data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

# Shape tells us the dimensions
print("Shape:", data.shape)  # Output: (10,)

# Size tells us total number of elements
print("Size:", data.size)    # Output: 10

# ndim tells us number of dimensions
print("Dimensions:", data.ndim)  # Output: 1

# dtype tells us the data type
print("Data type:", data.dtype)  # Output: int64

# Create an array with 36 months of sales data (1 through 36)
monthly_sales = np.arange(1, 37)  # 1, 2, 3, ..., 36
print("Original shape:", monthly_sales.shape)  # (36,)
print("First 10 months:", monthly_sales[:10])  # [1 2 3 4 5 6 7 8 9 10]

# Reshape to 3 years × 12 months
# We want 3 rows (years) and 12 columns (months)
yearly_sales = monthly_sales.reshape(3, 12)
print("\nReshaped array shape:", yearly_sales.shape)  # (3, 12)
print("Yearly sales array:")
print(yearly_sales)

Original shape: (36,)
First 10 months: [1 2 3 4 5 6 7 8 9 10]

Reshaped array shape: (3, 12)
Yearly sales array:
[[ 1  2  3  4  5  6  7  8  9 10 11 12]
 [13 14 15 16 17 18 19 20 21 22 23 24]
 [25 26 27 28 29 30 31 32 33 34 35 36]]

# Sample data: Daily temperatures for a week (°C)
temperatures = np.array([22, 24, 19, 25, 23, 21, 20])

# Calculate the mean
mean_temp = np.mean(temperatures)
print(f"Mean temperature: {mean_temp:.2f}°C")  # :.2f formats to 2 decimal places

# Let's verify by calculating manually
manual_mean = sum(temperatures) / len(temperatures)
print(f"Manual calculation: {manual_mean:.2f}°C")

Mean temperature: 22.00°C
Manual calculation: 22.00°C

# Add an outlier (a heat wave day!)
temperatures_with_outlier = np.array([22, 24, 19, 25, 23, 21, 20, 40])
print(f"Original mean: {np.mean(temperatures):.2f}°C")
print(f"Mean with outlier: {np.mean(temperatures_with_outlier):.2f}°C")
print(f"Change: {np.mean(temperatures_with_outlier) - np.mean(temperatures):.2f}°C")

Original mean: 22.00°C
Mean with outlier: 24.25°C
Change: 2.25°C

# Example 1: Odd number of values
temps_odd = np.array([19, 20, 21, 22, 23, 24, 25])
median_odd = np.median(temps_odd)
print(f"Dataset (odd): {temps_odd}")
print(f"Median (odd count): {median_odd}")  # Output: 22

# Example 2: Even number of values  
temps_even = np.array([19, 20, 21, 22, 23, 24])
median_even = np.median(temps_even)
print(f"\nDataset (even): {temps_even}")
print(f"Median (even count): {median_even}")  # Output: 21.5

# Same outlier example as before
temps_with_outlier = np.array([19, 20, 21, 22, 23, 24, 25, 40])  # Added 40 at the end
original_median = np.median([19, 20, 21, 22, 23, 24, 25])
new_median = np.median(temps_with_outlier)

print(f"Original median: {original_median}")
print(f"Median with outlier: {new_median}")
print(f"Change: {new_median - original_median}")

Original median: 22.0
Median with outlier: 22.5
Change: 0.5

# Let's use SciPy for mode (install with: pip install scipy)
from scipy import stats

# Survey responses on a scale of 1-5
survey_responses = np.array([1, 2, 2, 3, 3, 3, 4, 4, 5])
print(f"Survey responses: {survey_responses}")

mode_result = stats.mode(survey_responses)
print(f"Mode: {mode_result.mode[0]}")
print(f"How many times it appears: {mode_result.count[0]}")

Survey responses: [1 2 2 3 3 3 4 4 5]
Mode: 3
How many times it appears: 3

# Student grades with different assignment weights
grades = np.array([85, 92, 78, 95, 88])
print(f"Grades: {grades}")

# Weights for each assignment (must sum to 1 or 100%)
weights = np.array([0.1, 0.2, 0.2, 0.3, 0.2])  # Homework, Quiz, Midterm, Project, Final
print(f"Weights: {weights}")
print(f"Sum of weights: {np.sum(weights)}")  # Should be 1.0

# Calculate weighted average
weighted_avg = np.average(grades, weights=weights)
print(f"\nWeighted average: {weighted_avg:.2f}")

# Compare with regular mean
regular_mean = np.mean(grades)
print(f"Regular mean: {regular_mean:.2f}")
print(f"Difference: {weighted_avg - regular_mean:.2f}")

Grades: [85 92 78 95 88]
Weights: [0.1 0.2 0.2 0.3 0.2]
Sum of weights: 1.0

Weighted average: 89.30
Regular mean: 87.60
Difference: 1.70

# Let's verify step by step
weighted_sum = np.sum(grades * weights)
print(f"Weighted sum: {weighted_sum}")  # 88.6
print(f"NumPy weighted average: {np.average(grades, weights=weights)}")  # 89.3

# The CORRECT understanding:
# When you use np.average with weights, it does:
# 1. Multiply each value by its weight
# 2. Sum those products
# 3. Divide by the sum of weights

weights = np.array([1, 2, 2, 3, 2])  # These don't sum to 1, they're "importance" weights
weighted_avg = np.average(grades, weights=weights)
print(f"\nCorrect weighted average: {weighted_avg:.2f}")

# Manual calculation:
# (85×1 + 92×2 + 78×2 + 95×3 + 88×2) / (1+2+2+3+2)
numerator = 85*1 + 92*2 + 78*2 + 95*3 + 88*2  # 85 + 184 + 156 + 285 + 176 = 886
denominator = 1 + 2 + 2 + 3 + 2  # 10
manual_weighted = numerator / denominator  # 886 / 10 = 88.6

print(f"Manual calculation: {manual_weighted:.2f}")

Correct weighted average: 88.60
Manual calculation: 88.60

# Salary data for two companies
salaries_company_a = np.array([45000, 48000, 52000, 55000, 60000])
salaries_company_b = np.array([45000, 46000, 47000, 48000, 120000])

print("Company A salaries:", salaries_company_a)
print("Company B salaries:", salaries_company_b)

# Calculate range using ptp() - "peak to peak"
range_a = np.ptp(salaries_company_a)
range_b = np.ptp(salaries_company_b)

print(f"\nCompany A range: ${range_a:,}")
print(f"Company B range: ${range_b:,}")

# Manual calculation for comparison
manual_range_a = np.max(salaries_company_a) - np.min(salaries_company_a)
manual_range_b = np.max(salaries_company_b) - np.min(salaries_company_b)

print(f"\nManual range A: ${manual_range_a:,} (same as ptp)")
print(f"Manual range B: ${manual_range_b:,} (same as ptp)")

Company A salaries: [45000 48000 52000 55000 60000]
Company B salaries: [45000 46000 47000 48000 120000]

Company A range: $15,000
Company B range: $75,000

Manual range A: $15,000 (same as ptp)
Manual range B: $75,000 (same as ptp)

# Let's calculate variance step-by-step to understand it
data = np.array([2, 4, 6, 8, 10])
print(f"Data: {data}")

# Step 1: Find the mean
mean = np.mean(data)
print(f"\nStep 1 - Mean: {mean}")

# Step 2: Find differences from mean
differences = data - mean
print(f"Step 2 - Differences from mean: {differences}")

# Step 3: Square the differences (removes negatives, emphasizes large deviations)
squared_diffs = differences ** 2
print(f"Step 3 - Squared differences: {squared_diffs}")

# Step 4: Average the squared differences
variance_manual = np.mean(squared_diffs)
print(f"Step 4 - Average of squared differences: {variance_manual}")

# Step 5: Compare with NumPy's built-in
variance_numpy = np.var(data)
print(f"\nNumPy variance: {variance_numpy}")
print(f"Match? {np.isclose(variance_manual, variance_numpy)}")

Data: [2 4 6 8 10]

Step 1 - Mean: 6.0
Step 2 - Differences from mean: [-4. -2.  0.  2.  4.]
Step 3 - Squared differences: [16.  4.  0.  4. 16.]
Step 4 - Average of squared differences: 8.0

NumPy variance: 8.0
Match? True

# Same data as before
data = np.array([2, 4, 6, 8, 10])

# Calculate standard deviation
std_dev = np.std(data)
print(f"Standard deviation: {std_dev:.2f}")

# Verify it's the square root of variance
variance = np.var(data)
print(f"Variance: {variance:.2f}")
print(f"Square root of variance: {np.sqrt(variance):.2f}")
print(f"Match? {np.isclose(std_dev, np.sqrt(variance))}")

# Real-world example: Test scores
test_scores = np.array([78, 85, 92, 65, 88, 72, 95, 81])
mean_score = np.mean(test_scores)
std_score = np.std(test_scores)

print(f"\nTest scores: {test_scores}")
print(f"Mean score: {mean_score:.1f}")
print(f"Standard deviation: {std_score:.1f} points")
print(f"Most scores are within {std_score:.1f} points of {mean_score:.1f}")

Standard deviation: 2.83
Variance: 8.00
Square root of variance: 2.83
Match? True

Test scores: [78 85 92 65 88 72 95 81]
Mean score: 82.0
Standard deviation: 9.8 points
Most scores are within 9.8 points of 82.0

# Generate normally distributed data
np.random.seed(42)  # For reproducible results
normal_data = np.random.normal(100, 15, 10000)  # Mean=100, Std=15, 10000 points

mean = np.mean(normal_data)
std = np.std(normal_data)

# Calculate percentages within 1, 2, 3 standard deviations
within_1_std = np.sum((normal_data >= mean - std) & (normal_data <= mean + std)) / 10000 * 100
within_2_std = np.sum((normal_data >= mean - 2*std) & (normal_data <= mean + 2*std)) / 10000 * 100
within_3_std = np.sum((normal_data >= mean - 3*std) & (normal_data <= mean + 3*std)) / 10000 * 100

print(f"Mean: {mean:.2f}, Std: {std:.2f}")
print(f"Within 1 std: {within_1_std:.1f}% (theoretical: 68.3%)")
print(f"Within 2 std: {within_2_std:.1f}% (theoretical: 95.4%)")
print(f"Within 3 std: {within_3_std:.1f}% (theoretical: 99.7%)")

Mean: 100.02, Std: 14.98
Within 1 std: 68.3% (theoretical: 68.3%)
Within 2 std: 95.5% (theoretical: 95.4%)
Within 3 std: 99.7% (theoretical: 99.7%)

data = np.array([1, 2, 3, 4, 5])

# Population statistics (when you have ALL the data)
pop_var = np.var(data)  # Default: ddof=0
pop_std = np.std(data)  # Default: ddof=0

# Sample statistics (when you have a SAMPLE of a larger population)
sample_var = np.var(data, ddof=1)  # ddof=1 for sample
sample_std = np.std(data, ddof=1)  # ddof=1 for sample

print(f"Data: {data}")
print(f"\nPopulation variance: {pop_var:.4f}")
print(f"Sample variance: {sample_var:.4f}")
print(f"Difference: {sample_var - pop_var:.4f}")

print(f"\nPopulation std: {pop_std:.4f}")
print(f"Sample std: {sample_std:.4f}")
print(f"Difference: {sample_std - pop_std:.4f}")

Data: [1 2 3 4 5]

Population variance: 2.0000
Sample variance: 2.5000
Difference: 0.5000

Population std: 1.4142
Sample std: 1.5811
Difference: 0.1669

# House prices in a neighborhood (in thousands)
house_prices = np.array([250, 275, 300, 320, 350, 380, 400, 450, 500, 2000])
print(f"House prices: ${house_prices}K")
print(f"Note: $2000K ($2M) is an outlier - maybe a mansion!")

# Calculate quartiles
q1 = np.percentile(house_prices, 25)  # 25th percentile
q3 = np.percentile(house_prices, 75)  # 75th percentile
iqr = q3 - q1

print(f"\nQ1 (25th percentile): ${q1:.0f}K")
print(f"Q3 (75th percentile): ${q3:.0f}K")
print(f"IQR: ${iqr:.0f}K (middle 50% of houses)")

# Outlier detection using IQR
# Common rule: Values outside Q1 - 1.5×IQR or Q3 + 1.5×IQR are outliers
lower_bound = q1 - 1.5 * iqr
upper_bound = q3 + 1.5 * iqr

print(f"\nOutlier bounds: [${lower_bound:.0f}K, ${upper_bound:.0f}K]")

# Find outliers
outliers = house_prices[(house_prices < lower_bound) | (house_prices > upper_bound)]
print(f"Outliers detected: {outliers}")

# Count how many outliers
print(f"Number of outliers: {len(outliers)} out of {len(house_prices)} houses")

House prices: $[ 250  275  300  320  350  380  400  450  500 2000]K
Note: $2000K ($2M) is an outlier - maybe a mansion!

Q1 (25th percentile): $300K
Q3 (75th percentile): $425K
IQR: $125K (middle 50% of houses)

Outlier bounds: [$112K, $613K]
Outliers detected: [2000]
Number of outliers: 1 out of 10 houses

# Generate exam scores for a large class
np.random.seed(123)  # For reproducible random numbers
exam_scores = np.random.normal(75, 10, 1000)  # 1000 scores, mean=75, std=10

print(f"Generated {len(exam_scores)} exam scores")
print(f"Mean: {np.mean(exam_scores):.1f}, Std: {np.std(exam_scores):.1f}")

# Common percentiles
percentile_values = [25, 50, 75, 90, 95]
percentiles = np.percentile(exam_scores, percentile_values)

print("\nCommon Percentiles:")
print("-" * 30)
for p_value, p_score in zip(percentile_values, percentiles):
    print(f"{p_value}th percentile: {p_score:.1f}")

# What score do you need to beat 80% of students?
score_80th = np.percentile(exam_scores, 80)
print(f"\nTo beat 80% of students, you need: {score_80th:.1f}")

# Let's verify what percentile means
# Count how many scores are below the 80th percentile score
scores_below_80th = np.sum(exam_scores < score_80th)
percentage_below = scores_below_80th / len(exam_scores) * 100
print(f"Actual percentage below {score_80th:.1f}: {percentage_below:.1f}% (should be ~80%)")

Generated 1000 exam scores
Mean: 74.7, Std: 10.1

Common Percentiles:
------------------------------
25th percentile: 68.0
50th percentile: 74.6
75th percentile: 81.4
90th percentile: 87.6
95th percentile: 91.0

To beat 80% of students, you need: 83.2
Actual percentage below 83.2: 79.9% (should be ~80%)

# Get deciles (10%, 20%, ..., 90%) all at once
decile_percentiles = np.arange(10, 100, 10)  # [10, 20, 30, ..., 90]
deciles = np.percentile(exam_scores, decile_percentiles)

print("Deciles (every 10th percentile):")
print("-" * 40)
for percentile, score in zip(decile_percentiles, deciles):
    print(f"{percentile}th percentile: {score:.2f}")

# Visual representation
print("\nScore Distribution by Decile:")
print("-" * 40)
for i in range(len(deciles)):
    lower = deciles[i] if i == 0 else deciles[i-1]
    upper = deciles[i]
    print(f"Bottom {decile_percentiles[i]}%: scores < {upper:.1f}")

Deciles (every 10th percentile):
----------------------------------------
10th percentile: 61.82
20th percentile: 66.14
30th percentile: 69.55
40th percentile: 72.37
50th percentile: 74.64
60th percentile: 77.12
70th percentile: 79.88
80th percentile: 83.24
90th percentile: 87.64

Score Distribution by Decile:
----------------------------------------
Bottom 10%: scores < 61.8
Bottom 20%: scores < 66.1
Bottom 30%: scores < 69.6
Bottom 40%: scores < 72.4
Bottom 50%: scores < 74.6
Bottom 60%: scores < 77.1
Bottom 70%: scores < 79.9
Bottom 80%: scores < 83.2
Bottom 90%: scores < 87.6

# Hours studied vs. Exam scores
hours_studied = np.array([2, 5, 7, 10, 12, 15, 18, 20])
exam_scores = np.array([55, 65, 70, 75, 80, 85, 88, 90])

print(f"Hours studied: {hours_studied}")
print(f"Exam scores: {exam_scores}")

# Calculate covariance matrix
# np.cov() returns a 2x2 matrix for two variables
cov_matrix = np.cov(hours_studied, exam_scores)
print(f"\nCovariance matrix:\n{cov_matrix}")

# The covariance between hours and scores is at position [0, 1] or [1, 0]
covariance = cov_matrix[0, 1]
print(f"\nCovariance between hours studied and exam scores: {covariance:.2f}")

# What does this number mean?
print("\nInterpreting covariance:")
if covariance > 0:
    print("  Positive covariance: As hours increase, scores tend to increase")
elif covariance < 0:
    print("  Negative covariance: As hours increase, scores tend to decrease")
else:
    print("  Zero covariance: No linear relationship")

# Let's calculate covariance manually to understand it
mean_hours = np.mean(hours_studied)
mean_scores = np.mean(exam_scores)

# Step 1: Deviations from mean
dev_hours = hours_studied - mean_hours
dev_scores = exam_scores - mean_scores

# Step 2: Product of deviations
products = dev_hours * dev_scores

# Step 3: Average of products (for sample, divide by n-1)
n = len(hours_studied)
cov_manual = np.sum(products) / (n - 1)

print(f"\nManual calculation:")
print(f"  Mean hours: {mean_hours:.2f}")
print(f"  Mean scores: {mean_scores:.2f}")
print(f"  Deviations hours: {dev_hours}")
print(f"  Deviations scores: {dev_scores}")
print(f"  Products: {products}")
print(f"  Sum of products: {np.sum(products):.2f}")
print(f"  Covariance: {np.sum(products)} / ({n} - 1) = {cov_manual:.2f}")

Hours studied: [ 2  5  7 10 12 15 18 20]
Exam scores: [55 65 70 75 80 85 88 90]

Covariance matrix:
[[ 39.92857143 112.5       ]
 [112.5        142.85714286]]

Covariance between hours studied and exam scores: 112.50

Interpreting covariance:
  Positive covariance: As hours increase, scores tend to increase

Manual calculation:
  Mean hours: 11.12
  Mean scores: 76.00
  Deviations hours: [-9.125 -6.125 -4.125 -1.125  0.875  3.875  6.875  8.875]
  Deviations scores: [-21. -11.  -6.  -1.   4.   9.  12.  14.]
  Products: [191.625  67.375  24.75    1.125   3.5    34.875  82.5   124.25]
  Sum of products: 787.875
  Covariance: 787.875 / (8 - 1) = 112.50

# Same data
hours_studied = np.array([2, 5, 7, 10, 12, 15, 18, 20])
exam_scores = np.array([55, 65, 70, 75, 80, 85, 88, 90])

# Calculate correlation matrix
corr_matrix = np.corrcoef(hours_studied, exam_scores)
print(f"Correlation matrix:\n{corr_matrix}")

# Extract the correlation coefficient
correlation = corr_matrix[0, 1]
print(f"\nCorrelation coefficient: {correlation:.3f}")

# Interpretation guide
print("\nCorrelation Interpretation Guide:")
print("-" * 40)

if abs(correlation) >= 0.9:
    strength = "very strong"
elif abs(correlation) >= 0.7:
    strength = "strong"
elif abs(correlation) >= 0.5:
    strength = "moderate"
elif abs(correlation) >= 0.3:
    strength = "weak"
else:
    strength = "very weak or no"

if correlation > 0:
    direction = "positive"
elif correlation < 0:
    direction = "negative"
else:
    direction = "no"

print(f"This shows a {strength} {direction} correlation.")

print(f"\nIn plain English:")
if correlation > 0.7:
    print("  Students who study more tend to score much higher")
elif correlation > 0.3:
    print("  There's some relationship between studying and scores")
elif correlation > -0.3:
    print("  Studying time doesn't strongly predict exam scores")
else:
    print("  Students who study more tend to score lower (unlikely here!)")

# How correlation is calculated from covariance
covariance = np.cov(hours_studied, exam_scores)[0, 1]
std_hours = np.std(hours_studied, ddof=1)  # Sample std
std_scores = np.std(exam_scores, ddof=1)   # Sample std

correlation_manual = covariance / (std_hours * std_scores)
print(f"\nManual calculation:")
print(f"  Covariance: {covariance:.3f}")
print(f"  Std of hours: {std_hours:.3f}")
print(f"  Std of scores: {std_scores:.3f}")
print(f"  Correlation = {covariance:.3f} / ({std_hours:.3f} × {std_scores:.3f})")
print(f"  Correlation = {covariance:.3f} / {std_hours * std_scores:.3f}")
print(f"  Correlation = {correlation_manual:.3f}")
print(f"  Matches NumPy? {np.isclose(correlation, correlation_manual)}")

Correlation matrix:
[[1.         0.94491118]
 [0.94491118 1.        ]]

Correlation coefficient: 0.945

Correlation Interpretation Guide:
----------------------------------------
This shows a very strong positive correlation.

In plain English:
  Students who study more tend to score much higher

Manual calculation:
  Covariance: 112.500
  Std of hours: 6.319
  Std of scores: 11.952
  Correlation = 112.500 / (6.319 × 11.952)
  Correlation = 112.500 / 75.523
  Correlation = 0.945
  Matches NumPy? True

import matplotlib.pyplot as plt

# Generate example data with different correlations
np.random.seed(42)
n = 100

# Strong positive correlation
x = np.random.randn(n)
y_positive = x + np.random.randn(n) * 0.3
corr_positive = np.corrcoef(x, y_positive)[0, 1]

# Strong negative correlation  
y_negative = -x + np.random.randn(n) * 0.3
corr_negative = np.corrcoef(x, y_negative)[0, 1]

# No correlation
y_none = np.random.randn(n)
corr_none = np.corrcoef(np.random.randn(n), y_none)[0, 1]

# Create visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, y, corr, title in zip(axes, 
                              [y_positive, y_negative, y_none],
                              [corr_positive, corr_negative, corr_none],
                              [f"Positive Correlation (r={corr_positive:.2f})",
                               f"Negative Correlation (r={corr_negative:.2f})",
                               f"No Correlation (r={corr_none:.2f})"]):
    ax.scatter(x, y, alpha=0.6)
    ax.set_title(title)
    ax.set_xlabel("Variable X")
    ax.set_ylabel("Variable Y")
    ax.grid(True, alpha=0.3)

    # Add correlation interpretation
    if abs(corr) > 0.7:
        strength = "Strong"
    elif abs(corr) > 0.3:
        strength = "Moderate"
    else:
        strength = "Weak/None"

    ax.text(0.05, 0.95, f"{strength} relationship", 
            transform=ax.transAxes, fontsize=10,
            verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.tight_layout()
plt.show()

# Real data often has missing values
# Let's say we have daily sales data, but some days data wasn't recorded
sales_data = np.array([1200, 1350, np.nan, 1420, np.nan, 1550, 1480])
print(f"Sales data with missing values: {sales_data}")
print(f"Type of np.nan: {type(np.nan)}")

# Problems with NaN in regular operations
print("\nProblems with NaN in calculations:")
print(f"Mean with regular np.mean(): {np.mean(sales_data)}")  # Returns nan!
print(f"Sum with regular np.sum(): {np.sum(sales_data)}")     # Also nan!
print(f"Max with regular np.max(): {np.max(sales_data)}")     # Even max returns nan!

# The issue: Any operation with NaN propagates NaN
# NaN is like a virus - it infects any calculation it touches

# Solution: Use nan-aware functions
print("\nSolution: Use nan-aware functions:")
print(f"Nan-safe mean: {np.nanmean(sales_data):.2f}")
print(f"Nan-safe sum: {np.nansum(sales_data):.0f}")
print(f"Nan-safe standard deviation: {np.nanstd(sales_data):.2f}")
print(f"Nan-safe max: {np.nanmax(sales_data)}")
print(f"Nan-safe min: {np.nanmin(sales_data)}")

# Count non-NaN values
valid_count = np.count_nonzero(~np.isnan(sales_data))
total_count = len(sales_data)
print(f"\nData completeness: {valid_count}/{total_count} values are valid ({valid_count/total_count*100:.1f}%)")

# Find which positions have NaN
nan_positions = np.where(np.isnan(sales_data))[0]
print(f"NaN at positions: {nan_positions}")

# Replace NaN with a value (like mean)
sales_filled = sales_data.copy()  # Make a copy to avoid changing original
mean_sales = np.nanmean(sales_data)
sales_filled[np.isnan(sales_filled)] = mean_sales
print(f"\nData with NaN filled with mean ({mean_sales:.0f}): {sales_filled.astype(int)}")

Sales data with missing values: [1200. 1350.   nan 1420.   nan 1550. 1480.]
Type of np.nan: <class 'float'>

Problems with NaN in calculations:
Mean with regular np.mean(): nan
Sum with regular np.sum(): nan
Max with regular np.max(): nan

Solution: Use nan-aware functions:
Nan-safe mean: 1400.00
Nan-safe sum: 7000
Nan-safe standard deviation: 130.38
Nan-safe max: 1550.0
Nan-safe min: 1200.0

Data completeness: 5/7 values are valid (71.4%)
NaN at positions: [2 4]

Data with NaN filled with mean (1400): [1200 1350 1400 1420 1400 1550 1480]

# Be careful with integer vs float arrays
int_array = np.array([1, 2, 3, 4, 5])
float_array = np.array([1., 2., 3., 4., 5.])

print(f"Integer array: {int_array}, dtype: {int_array.dtype}")
print(f"Float array: {float_array}, dtype: {float_array.dtype}")

# Mean of integers
int_mean = np.mean(int_array)
print(f"\nMean of ints: {int_mean} - type: {int_mean.dtype}")

# Division with integers can be surprising
print(f"\nInt division by 2: {int_array / 2}")  # Becomes float!
print(f"Floor division: {int_array // 2}")      # Keeps as int (truncates)

# What happens with different operations?
print(f"\nInteger operations:")
print(f"  int + int: {(int_array + 2).dtype}")      # Stays int
print(f"  int * float: {(int_array * 2.0).dtype}")  # Becomes float
print(f"  int / int: {(int_array / int_array).dtype}")  # Becomes float!

# Convert between types
float_version = int_array.astype(float)
int_version = float_array.astype(int)  # Truncates (1.9 becomes 1!)
print(f"\nType conversion:")
print(f"  int to float: {float_version}, dtype: {float_version.dtype}")
print(f"  float to int: {int_version}, dtype: {int_version.dtype}")

# Be careful with large numbers!
large_int = np.array([9999999999], dtype=np.int32)  # 32-bit int
print(f"\nLarge int (32-bit): {large_int}")
print(f"  Might overflow! Max 32-bit int: {2**31 - 1:,}")

# Safe conversion
safe_float = large_int.astype(np.float64)  # 64-bit float can handle it
print(f"  Safe as float64: {safe_float}")

Integer array: [1 2 3 4 5], dtype: int64
Float array: [1. 2. 3. 4. 5.], dtype: float64

Mean of ints: 3.0 - type: float64

Int division by 2: [0.5 1.  1.5 2.  2.5]
Floor division: [0 1 1 2 2]

Integer operations:
  int + int: int64
  int * float: float64
  int / int: float64

Type conversion:
  int to float: [1. 2. 3. 4. 5.], dtype: float64
  float to int: [1 2 3 4 5], dtype: int64

Large int (32-bit): [1410065407]
  Might overflow! Max 32-bit int: 2,147,483,647
  Safe as float64: [1410065407.]

# For huge datasets, memory and speed matter
# Let's create a large dataset to demonstrate
large_data = np.random.randn(10_000_000)  # 10 million random numbers!
print(f"Created array with {large_data.size:,} elements")
print(f"Memory usage: {large_data.nbytes / 1_000_000:.1f} MB")
print(f"Data type: {large_data.dtype}")
print(f"Shape: {large_data.shape}")

# Time different ways to calculate mean
import time

print("\nTiming different mean calculations:")
print("-" * 40)

# Method 1: Python's built-in sum (slow)
start = time.time()
mean_slow = sum(large_data) / len(large_data)  # Python's sum
time_slow = time.time() - start
print(f"Python sum time: {time_slow:.3f} seconds")

# Method 2: NumPy's mean (fast)
start = time.time()
mean_fast = np.mean(large_data)  # NumPy's vectorized mean
time_fast = time.time() - start
print(f"NumPy mean time: {time_fast:.3f} seconds")

print(f"\nNumPy is {time_slow/time_fast:.1f}x faster!")
print(f"Results match? {np.isclose(mean_slow, mean_fast)}")

# Memory optimization: Use smaller data types
print("\n" + "="*50)
print("Memory Optimization Examples:")
print("="*50)

# Float64 (default)
float64_data = np.ones(1_000_000, dtype=np.float64)
print(f"\nFloat64 array (1M elements):")
print(f"  Memory: {float64_data.nbytes / 1_000_000:.1f} MB")
print(f"  Precision: ~15 decimal digits")

# Float32 (half the memory)
float32_data = np.ones(1_000_000, dtype=np.float32)
print(f"\nFloat32 array (1M elements):")
print(f"  Memory: {float32_data.nbytes / 1_000_000:.1f} MB")
print(f"  Precision: ~7 decimal digits")
print(f"  Memory saved: {(float64_data.nbytes - float32_data.nbytes) / 1_000_000:.1f} MB")

# Int8 (even smaller)
int8_data = np.ones(1_000_000, dtype=np.int8)
print(f"\nInt8 array (1M elements):")
print(f"  Memory: {int8_data.nbytes / 1_000_000:.1f} MB")
print(f"  Range: -128 to 127")

# When to use which?
print("\n" + "="*50)
print("Guidelines for choosing data types:")
print("="*50)
print("1. Use float64 for precise calculations (default)")
print("2. Use float32 for large datasets where precision less critical")
print("3. Use integers when dealing with counts or whole numbers")
print("4. Use smallest type that fits your data range")
print("5. Check memory with `.nbytes` property")

Created array with 10,000,000 elements
Memory usage: 80.0 MB
Data type: float64
Shape: (10000000,)

Timing different mean calculations:
----------------------------------------
Python sum time: 0.423 seconds
NumPy mean time: 0.008 seconds

NumPy is 52.9x faster!
Results match? True

==================================================
Memory Optimization Examples:
==================================================

Float64 array (1M elements):
  Memory: 8.0 MB
  Precision: ~15 decimal digits

Float32 array (1M elements):
  Memory: 4.0 MB
  Precision: ~7 decimal digits
  Memory saved: 4.0 MB

Int8 array (1M elements):
  Memory: 1.0 MB
  Range: -128 to 127

==================================================
Guidelines for choosing data types:
==================================================
1. Use float64 for precise calculations (default)
2. Use float32 for large datasets where precision less critical
3. Use integers when dealing with counts or whole numbers
4. Use smallest type that fits your data range
5. Check memory with `.nbytes` property

# Monthly sales for 3 products across 4 quarters
# Rows = products, Columns = quarters
sales_data = np.array([
    [120, 150, 130, 140],  # Product A quarterly sales
    [80, 90, 85, 95],      # Product B quarterly sales  
    [200, 220, 210, 230]   # Product C quarterly sales
])

print("Sales Data (3 products × 4 quarters):")
print(sales_data)
print(f"\nShape: {sales_data.shape}")  # (3 products, 4 quarters)
print(f"Size: {sales_data.size}")      # 3 × 4 = 12 elements
print(f"Dimensions: {sales_data.ndim}") # 2 dimensions

# Mean of each product (across quarters) - axis=1
print("\n1. Average quarterly sales PER PRODUCT (axis=1):")
print("-" * 50)
product_means = np.mean(sales_data, axis=1)
for i, mean in enumerate(product_means, 1):
    print(f"  Product {i}: ${mean:.1f} average per quarter")

# How axis=1 works:
# For each row (product), calculate mean across columns (quarters)
# Row 1 (Product A): mean(120, 150, 130, 140) = 135
# Row 2 (Product B): mean(80, 90, 85, 95) = 87.5
# Row 3 (Product C): mean(200, 220, 210, 230) = 215

# Mean of each quarter (across products) - axis=0
print("\n2. Average sales PER QUARTER (axis=0):")
print("-" * 50)
quarter_means = np.mean(sales_data, axis=0)
quarters = ['Q1', 'Q2', 'Q3', 'Q4']
for q, mean in zip(quarters, quarter_means):
    print(f"  {q}: ${mean:.1f} average across all products")

# How axis=0 works:
# For each column (quarter), calculate mean across rows (products)
# Column 1 (Q1): mean(120, 80, 200) = 133.33
# Column 2 (Q2): mean(150, 90, 220) = 153.33
# Column 3 (Q3): mean(130, 85, 210) = 141.67
# Column 4 (Q4): mean(140, 95, 230) = 155

# Overall mean (all elements)
print("\n3. OVERALL average (axis=None):")
print("-" * 50)
overall_mean = np.mean(sales_data)  # axis=None is default
print(f"  Overall: ${overall_mean:.2f}")

# How overall mean works:
# mean(120, 150, 130, 140, 80, 90, 85, 95, 200, 220, 210, 230) = 145.83

# Variance along different axes
print("\n4. Variance analysis:")
print("-" * 50)
print(f"Variance per product (across quarters): {np.var(sales_data, axis=1)}")
print(f"Variance per quarter (across products): {np.var(sales_data, axis=0)}")
print(f"Overall variance: {np.var(sales_data):.2f}")

# Understanding the axis parameter visually
print("\n" + "="*60)
print("VISUAL GUIDE to the axis parameter:")
print("="*60)
print("Array: 3 rows × 4 columns")
print("[[120, 150, 130, 140]")
print(" [ 80,  90,  85,  95]")
print(" [200, 220, 210, 230]]")
print("\naxis=0 → Apply vertically (column-wise)")
print("  Think: 'Collapse rows, keep columns'")
print("  Result: 4 values (one per column)")
print("\naxis=1 → Apply horizontally (row-wise)")
print("  Think: 'Collapse columns, keep rows'")
print("  Result: 3 values (one per row)")
print("\naxis=None → Apply to all elements")
print("  Think: 'Collapse everything'")
print("  Result: 1 value")

Sales Data (3 products × 4 quarters):
[[120 150 130 140]
 [ 80  90  85  95]
 [200 220 210 230]]

Shape: (3, 4)
Size: 12
Dimensions: 2

1. Average quarterly sales PER PRODUCT (axis=1):
--------------------------------------------------
  Product 1: $135.0 average per quarter
  Product 2: $87.5 average per quarter
  Product 3: $215.0 average per quarter

2. Average sales PER QUARTER (axis=0):
--------------------------------------------------
  Q1: $133.3 average across all products
  Q2: $153.3 average across all products
  Q3: $141.7 average across all products
  Q4: $155.0 average across all products

3. OVERALL average (axis=None):
--------------------------------------------------
  Overall: $145.83

4. Variance analysis:
--------------------------------------------------
Variance per product (across quarters): [116.66666667  29.16666667 116.66666667]
Variance per quarter (across products): [2600.         2600.         2600.         2600.        ]
Overall variance: 2600.00

============================================================
VISUAL GUIDE to the axis parameter:
============================================================
Array: 3 rows × 4 columns
[[120, 150, 130, 140]
 [ 80,  90,  85,  95]
 [200, 220, 210, 230]]

axis=0 → Apply vertically (column-wise)
  Think: 'Collapse rows, keep columns'
  Result: 4 values (one per column)

axis=1 → Apply horizontally (row-wise)
  Think: 'Collapse columns, keep rows'
  Result: 3 values (one per row)

axis=None → Apply to all elements
  Think: 'Collapse everything'
  Result: 1 value

# Three business metrics over 5 months
sales = np.array([100, 120, 130, 150, 170])        # in thousands
marketing = np.array([50, 55, 60, 65, 70])         # in thousands
website_traffic = np.array([1000, 1100, 1200, 1300, 1400])  # visitors

print("Monthly Business Data:")
print("Month  Sales($K)  Marketing($K)  Traffic")
print("-" * 40)
for i in range(len(sales)):
    print(f"{i+1:2d}     {sales[i]:5d}        {marketing[i]:5d}       {website_traffic[i]:5d}")

# Combine into a 2D array (5 months × 3 variables)
# Each column is a variable, each row is a month
data_matrix = np.column_stack([sales, marketing, website_traffic])
print(f"\nData matrix shape: {data_matrix.shape}")  # (5, 3)

# Calculate correlation matrix
# rowvar=False means columns are variables (rows are observations)
corr_matrix = np.corrcoef(data_matrix, rowvar=False)

print("\nCorrelation Matrix:")
print("=" * 50)
print("               Sales   Marketing   Traffic")
print("-" * 50)
labels = ['Sales', 'Marketing', 'Traffic']
for i in range(3):
    print(f"{labels[i]:10s}", end="")
    for j in range(3):
        print(f"   {corr_matrix[i, j]:7.3f} ", end="")
    print()

print("\n" + "=" * 50)
print("INTERPRETATION:")
print("=" * 50)

# Sales vs Marketing
corr_sales_marketing = corr_matrix[0, 1]
print(f"\n1. Sales vs Marketing: r = {corr_sales_marketing:.3f}")
if corr_sales_marketing > 0.7:
    print("   Strong positive correlation")
    print("   When marketing spend increases, sales tend to increase")
elif corr_sales_marketing > 0:
    print("   Positive correlation")
    print("   Marketing and sales move together")

# Sales vs Traffic
corr_sales_traffic = corr_matrix[0, 2]
print(f"\n2. Sales vs Website Traffic: r = {corr_sales_traffic:.3f}")
if corr_sales_traffic > 0.9:
    print("   Very strong positive correlation")
    print("   Website traffic is an excellent predictor of sales")
elif corr_sales_traffic > 0.7:
    print("   Strong correlation")
    print("   More traffic generally means more sales")

# Marketing vs Traffic
corr_marketing_traffic = corr_matrix[1, 2]
print(f"\n3. Marketing vs Website Traffic: r = {corr_marketing_traffic:.3f}")
if corr_marketing_traffic > 0.7:
    print("   Strong correlation")
    print("   Marketing efforts drive website traffic")

# Diagonal is always 1
print(f"\n4. Diagonal values are all 1.0")
print("   A variable perfectly correlates with itself (obviously!)")

# Matrix is symmetric
print(f"\n5. Matrix is symmetric")
print("   corr(Sales, Marketing) = corr(Marketing, Sales) = {corr_matrix[0, 1]:.3f}")
print(f"   Upper triangle = Lower triangle")

Monthly Business Data:
Month  Sales($K)  Marketing($K)  Traffic
----------------------------------------
 1       100           50       1000
 2       120           55       1100
 3       130           60       1200
 4       150           65       1300
 5       170           70       1400

Data matrix shape: (5, 3)

Correlation Matrix:
==================================================
               Sales   Marketing   Traffic
--------------------------------------------------
Sales          1.000     1.000     1.000 
Marketing      1.000     1.000     1.000 
Traffic        1.000     1.000     1.000 

==================================================
INTERPRETATION:
==================================================

1. Sales vs Marketing: r = 1.000
   Strong positive correlation
   When marketing spend increases, sales tend to increase

2. Sales vs Website Traffic: r = 1.000
   Very strong positive correlation
   Website traffic is an excellent predictor of sales

3. Marketing vs Website Traffic: r = 1.000
   Strong correlation
   Marketing efforts drive website traffic

4. Diagonal values are all 1.0
   A variable perfectly correlates with itself (obviously!)

5. Matrix is symmetric
   corr(Sales, Marketing) = corr(Marketing, Sales) = 1.000
   Upper triangle = Lower triangle

# More realistic data with some randomness
np.random.seed(42)
sales = np.array([100, 120, 130, 150, 170]) + np.random.randn(5) * 10
marketing = np.array([50, 55, 60, 65, 70]) + np.random.randn(5) * 5
website_traffic = np.array([1000, 1100, 1200, 1300, 1400]) + np.random.randn(5) * 50

data_matrix = np.column_stack([sales, marketing, website_traffic])
corr_matrix = np.corrcoef(data_matrix, rowvar=False)

print("\nMore Realistic Correlation Matrix (with random noise):")
print("=" * 60)
print("               Sales   Marketing   Traffic")
print("-" * 60)
labels = ['Sales', 'Marketing', 'Traffic']
for i in range(3):
    print(f"{labels[i]:10s}", end="")
    for j in range(3):
        print(f"   {corr_matrix[i, j]:7.3f} ", end="")
    print()

print("\nInterpretation of correlation coefficients:")
print(f"- Sales vs Marketing: {corr_matrix[0, 1]:.3f} → Very strong positive")
print(f"- Sales vs Traffic: {corr_matrix[0, 2]:.3f} → Strong positive")  
print(f"- Marketing vs Traffic: {corr_matrix[1, 2]:.3f} → Moderate positive")

# Calculate R-squared (coefficient of determination)
print("\nR-squared (explained variance):")
print(f"- Marketing explains {corr_matrix[0, 1]**2 * 100:.1f}% of sales variance")
print(f"- Traffic explains {corr_matrix[0, 2]**2 * 100:.1f}% of sales variance")

More Realistic Correlation Matrix (with random noise):
============================================================
               Sales   Marketing   Traffic
------------------------------------------------------------
Sales          1.000     0.976     0.958 
Marketing      0.976     1.000     0.866 
Traffic        0.958     0.866     1.000 

Interpretation of correlation coefficients:
- Sales vs Marketing: 0.976 → Very strong positive
- Sales vs Traffic: 0.958 → Strong positive
- Marketing vs Traffic: 0.866 → Moderate positive

R-squared (explained variance):
- Marketing explains 95.3% of sales variance
- Traffic explains 91.8% of sales variance

# First, install SciPy if you haven't: pip install scipy
from scipy import stats

print("Skewness and Kurtosis Examples")
print("=" * 60)

# Generate different types of distributions
np.random.seed(123)
normal_data = np.random.normal(0, 1, 1000)  # Normal distribution (bell curve)
skewed_right = np.random.exponential(1, 1000)  # Right-skewed (like income)
uniform_data = np.random.uniform(-1, 1, 1000)  # Uniform distribution (flat)

print("\n1. Normal Distribution (Bell Curve):")
print(f"   Mean: {np.mean(normal_data):.3f}")
print(f"   Median: {np.median(normal_data):.3f}")
print(f"   Skewness: {stats.skew(normal_data):.3f} (should be near 0)")
print(f"   Kurtosis: {stats.kurtosis(normal_data):.3f} (should be near 0)")

print("\n2. Right-Skewed Distribution (like income):")
print(f"   Mean: {np.mean(skewed_right):.3f}")
print(f"   Median: {np.median(skewed_right):.3f}")
print(f"   Skewness: {stats.skew(skewed_right):.3f} (positive = right skew)")
print(f"   Kurtosis: {stats.kurtosis(skewed_right):.3f}")

print("\n3. Uniform Distribution (flat):")
print(f"   Mean: {np.mean(uniform_data):.3f}")
print(f"   Median: {np.median(uniform_data):.3f}")
print(f"   Skewness: {stats.skew(uniform_data):.3f} (near 0)")
print(f"   Kurtosis: {stats.kurtosis(uniform_data):.3f} (negative = flatter than normal)")

# What skewness tells us
print("\n" + "="*60)
print("INTERPRETING SKEWNESS:")
print("="*60)
print("Skewness measures asymmetry:")
print("• 0: Symmetric (normal distribution)")
print("• Positive: Right-skewed (tail on right)")
print("  - Mean > Median")
print("  - Example: Income (few very high incomes)")
print("• Negative: Left-skewed (tail on left)")
print("  - Mean < Median")
print("  - Example: Exam scores (most do well, few fail)")

# What kurtosis tells us  
print("\n" + "="*60)
print("INTERPRETING KURTOSIS:")
print("="*60)
print("Kurtosis measures 'tailedness' (scipy gives excess kurtosis):")
print("• 0: Normal tails (mesokurtic)")
print("• Positive: Heavy tails (leptokurtic)")
print("  - More outliers than normal")
print("  - Example: Financial returns (crashes and booms)")
print("• Negative: Light tails (platykurtic)")
print("  - Fewer outliers than normal")
print("  - Example: Uniform distribution")

# Real example: Income data is usually right-skewed
print("\n" + "="*60)
print("REAL-WORLD EXAMPLE: Income Distribution")
print("="*60)

# Simulate income data (right-skewed)
# Most people earn moderate incomes, few earn very high incomes
incomes = np.random.lognormal(mean=10.5, sigma=0.8, size=1000)
incomes = np.round(incomes, -2)  # Round to nearest $100

print(f"Sample size: {len(incomes)} people")
print(f"Mean income: ${np.mean(incomes):,.0f}")
print(f"Median income: ${np.median(incomes):,.0f}")
print(f"Difference (mean - median): ${np.mean(incomes) - np.median(incomes):,.0f}")
print(f"Skewness: {stats.skew(incomes):.3f} (positive = right-skewed)")
print(f"Kurtosis: {stats.kurtosis(incomes):.3f}")

if stats.skew(incomes) > 0.5:
    print("\nConclusion: Income is right-skewed")
    print("- Most people earn near the median")
    print("- A few high earners pull the mean upward")
    print("- Median is better measure of 'typical' income")

Skewness and Kurtosis Examples
============================================================

1. Normal Distribution (Bell Curve):
   Mean: -0.011
   Median: -0.010
   Skewness: 0.036 (should be near 0)
   Kurtosis: 0.039 (should be near 0)

2. Right-Skewed Distribution (like income):
   Mean: 0.984
   Median: 0.669
   Skewness: 1.957 (positive = right skew)
   Kurtosis: 5.426

3. Uniform Distribution (flat):
   Mean: -0.002
   Median: 0.005
   Skewness: -0.010 (near 0)
   Kurtosis: -1.195 (negative = flatter than normal)

============================================================
INTERPRETING SKEWNESS:
============================================================
Skewness measures asymmetry:
* 0: Symmetric (normal distribution)
* Positive: Right-skewed (tail on right)
  - Mean > Median
  - Example: Income (few very high incomes)
* Negative: Left-skewed (tail on left)
  - Mean < Median
  - Example: Exam scores (most do well, few fail)

============================================================
INTERPRETING KURTOSIS:
============================================================
Kurtosis measures 'tailedness' (scipy gives excess kurtosis):
* 0: Normal tails (mesokurtic)
* Positive: Heavy tails (leptokurtic)
  - More outliers than normal
  - Example: Financial returns (crashes and booms)
* Negative: Light tails (platykurtic)
  - Fewer outliers than normal
  - Example: Uniform distribution

============================================================
REAL-WORLD EXAMPLE: Income Distribution
============================================================
Sample size: 1000 people
Mean income: $38,717
Median income: $36,100
Difference (mean - median): $2,617
Skewness: 1.053 (positive = right-skewed)
Kurtosis: 1.870

Conclusion: Income is right-skewed
- Most people earn near the median
- A few high earners pull the mean upward
- Median is better measure of 'typical' income

# Running sales total over days
daily_sales = np.array([100, 150, 200, 180, 220])
days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri']

print("Daily Sales Data:")
print("-" * 30)
for day, sales in zip(days, daily_sales):
    print(f"{day}: ${sales}")

# Cumulative sum
cumulative_sales = np.cumsum(daily_sales)
print("\nCumulative Sales:")
print("-" * 30)
for day, cum_sales in zip(days, cumulative_sales):
    print(f"{day}: ${cum_sales}")

# Running average
running_avg = np.cumsum(daily_sales) / np.arange(1, len(daily_sales) + 1)
print("\nRunning Average:")
print("-" * 30)
for day, avg in zip(days, running_avg):
    print(f"{day}: ${avg:.2f}")

# Cumulative product (for growth rates)
print("\n" + "="*50)
print("Example: Investment Growth")
print("="*50)

# Daily returns (as multipliers, e.g., 1.01 = 1% gain)
daily_returns = np.array([1.01, 0.99, 1.02, 1.005, 1.015])
cumulative_return = np.cumprod(daily_returns)

print("Day  Daily Return  Cumulative Return  Total Return")
print("-" * 50)
for i in range(len(daily_returns)):
    total_return = (cumulative_return[i] - 1) * 100
    print(f"{i+1:2d}     {daily_returns[i]:6.3f}          {cumulative_return[i]:6.3f}         {total_return:6.2f}%")

# Cumulative maximum (tracking record highs)
print("\n" + "="*50)
print("Example: Tracking Record High Temperatures")
print("="*50)

weekly_temps = np.array([72, 75, 70, 78, 76, 80, 79])
record_highs = np.maximum.accumulate(weekly_temps)

print("Day  Temp  Record High")
print("-" * 25)
days_of_week = ['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat']
for day, temp, record in zip(days_of_week, weekly_temps, record_highs):
    if temp == record:
        marker = "← NEW RECORD!"
    else:
        marker = ""
    print(f"{day:3s}  {temp:3d}°F     {record:3d}°F  {marker}")

Daily Sales Data:
------------------------------
Mon: $100
Tue: $150
Wed: $200
Thu: $180
Fri: $220

Cumulative Sales:
------------------------------
Mon: $100
Tue: $250
Wed: $450
Thu: $630
Fri: $850

Running Average:
------------------------------
Mon: $100.00
Tue: $125.00
Wed: $150.00
Thu: $157.50
Fri: $170.00

==================================================
Example: Investment Growth
==================================================
Day  Daily Return  Cumulative Return  Total Return
--------------------------------------------------
 1     1.010          1.010          1.00%
 2     0.990          1.000          0.00%
 3     1.020          1.020          2.00%
 4     1.005          1.025          2.52%
 5     1.015          1.040          4.02%

==================================================
Example: Tracking Record High Temperatures
==================================================
Day  Temp  Record High
-------------------------
Sun   72°F     72°F  ← NEW RECORD!
Mon   75°F     75°F  ← NEW RECORD!
Tue   70°F     75°F  
Wed   78°F     78°F  ← NEW RECORD!
Thu   76°F     78°F  
Fri   80°F     80°F  ← NEW RECORD!
Sat   79°F     80°F

def data_summary(data, name="Data", round_digits=4):
    """
    Print comprehensive statistical summary of data.

    Parameters:
    -----------
    data : numpy array
        The data to analyze
    name : str
        Name to display for the data
    round_digits : int
        Number of decimal places to round to
    """
    print(f"\n{'='*60}")
    print(f"COMPREHENSIVE STATISTICAL SUMMARY: {name}")
    print(f"{'='*60}")

    # Handle NaN values
    has_nan = np.any(np.isnan(data))
    if has_nan:
        data_clean = data[~np.isnan(data)]
        print(f"⚠️  Data contains {np.sum(np.isnan(data))} NaN values")
        print(f"   Analysis uses {len(data_clean)} non-NaN values")
    else:
        data_clean = data

    # Basic counts
    print(f"\n📊 BASIC INFORMATION:")
    print(f"   Count: {len(data_clean):,}")
    print(f"   NaN values: {np.sum(np.isnan(data))}")
    print(f"   Data type: {data.dtype}")
    print(f"   Shape: {data.shape}")

    # Measures of central tendency
    print(f"\n🎯 CENTRAL TENDENCY:")
    print(f"   Mean: {np.nanmean(data):.{round_digits}f}")
    print(f"   Median: {np.nanmedian(data):.{round_digits}f}")

    # Try to calculate mode (requires SciPy)
    try:
        from scipy import stats
        mode_result = stats.mode(data_clean, keepdims=True)
        if len(mode_result.mode) > 0:
            print(f"   Mode: {mode_result.mode[0]:.{round_digits}f}")
            print(f"   Mode count: {mode_result.count[0]} occurrences")
    except ImportError:
        print("   Mode: (Install SciPy for mode calculation)")

    # Measures of spread
    print(f"\n📈 SPREAD/VARIABILITY:")
    print(f"   Standard Deviation: {np.nanstd(data):.{round_digits}f}")
    print(f"   Variance: {np.nanvar(data):.{round_digits}f}")
    print(f"   Range: {np.nanmax(data) - np.nanmin(data):.{round_digits}f}")
    print(f"   IQR: {np.nanpercentile(data, 75) - np.nanpercentile(data, 25):.{round_digits}f}")

    # Minimum and maximum
    print(f"\n📉 MIN/MAX VALUES:")
    print(f"   Minimum: {np.nanmin(data):.{round_digits}f}")
    print(f"   25th Percentile (Q1): {np.nanpercentile(data, 25):.{round_digits}f}")
    print(f"   50th Percentile (Median): {np.nanpercentile(data, 50):.{round_digits}f}")
    print(f"   75th Percentile (Q3): {np.nanpercentile(data, 75):.{round_digits}f}")
    print(f"   Maximum: {np.nanmax(data):.{round_digits}f}")

    # Skewness and kurtosis (if SciPy available)
    try:
        from scipy import stats
        skew_val = stats.skew(data_clean)
        kurt_val = stats.kurtosis(data_clean)

        print(f"\n📊 DISTRIBUTION SHAPE:")
        print(f"   Skewness: {skew_val:.{round_digits}f}")

        # Interpret skewness
        if abs(skew_val) < 0.5:
            skew_text = "Approximately symmetric"
        elif skew_val > 0:
            skew_text = "Right-skewed (tail on right)"
        else:
            skew_text = "Left-skewed (tail on left)"
        print(f"     → {skew_text}")

        print(f"   Kurtosis (excess): {kurt_val:.{round_digits}f}")

        # Interpret kurtosis
        if abs(kurt_val) < 0.5:
            kurt_text = "Similar to normal distribution"
        elif kurt_val > 0:
            kurt_text = "Heavy tails (more outliers than normal)"
        else:
            kurt_text = "Light tails (fewer outliers than normal)"
        print(f"     → {kurt_text}")

    except ImportError:
        print(f"\n📊 DISTRIBUTION SHAPE:")
        print("   Install SciPy for skewness and kurtosis")

    # Outlier detection using IQR method
    q1 = np.nanpercentile(data, 25)
    q3 = np.nanpercentile(data, 75)
    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr

    outliers = data_clean[(data_clean < lower_bound) | (data_clean > upper_bound)]

    print(f"\n🔍 OUTLIER DETECTION (IQR method):")
    print(f"   Lower bound: {lower_bound:.{round_digits}f}")
    print(f"   Upper bound: {upper_bound:.{round_digits}f}")
    print(f"   Potential outliers: {len(outliers)} values")
    if len(outliers) > 0:
        print(f"   Outlier values: {outliers}")

    # Data quality checks
    print(f"\n✅ DATA QUALITY CHECKS:")

    # Check for constant data
    if np.nanstd(data) == 0:
        print("   ⚠️  Warning: Zero standard deviation (constant data)")

    # Check for infinite values
    if np.any(np.isinf(data)):
        print("   ⚠️  Warning: Data contains infinite values")

    # Check range
    data_range = np.nanmax(data) - np.nanmin(data)
    if data_range == 0:
        print("   ⚠️  Warning: Data range is zero")

    print(f"{'='*60}")

# Test the function
print("TESTING THE COMPREHENSIVE SUMMARY FUNCTION")
print("=" * 60)

# Test with normal data
normal_test = np.random.normal(100, 15, 1000)
data_summary(normal_test, "Normal Distribution Test", round_digits=3)

# Test with real-world-like data (with NaN)
print("\n\n" + "="*60)
print("TEST WITH REAL-WORLD DATA (WITH NaN VALUES)")
print("=" * 60)

real_world_data = np.array([10, 12, 15, np.nan, 18, 20, 22, 100, np.nan, 25])
data_summary(real_world_data, "Real World Data with Outliers", round_digits=2)

TESTING THE COMPREHENSIVE SUMMARY FUNCTION
============================================================

============================================================
COMPREHENSIVE STATISTICAL SUMMARY: Normal Distribution Test
============================================================

📊 BASIC INFORMATION:
   Count: 1,000
   NaN values: 0
   Data type: float64
   Shape: (1000,)

🎯 CENTRAL TENDENCY:
   Mean: 99.799
   Median: 99.756
   Mode: 75.228
   Mode count: 3 occurrences

📈 SPREAD/VARIABILITY:
   Standard Deviation: 15.024
   Variance: 225.711
   Range: 87.564
   IQR: 20.388

📉 MIN/MAX VALUES:
   Minimum: 52.434
   25th Percentile (Q1): 89.790
   50th Percentile (Median): 99.756
   75th Percentile (Q3): 110.178
   Maximum: 139.998

📊 DISTRIBUTION SHAPE:
   Skewness: 0.004
     → Approximately symmetric
   Kurtosis (excess): -0.140
     → Similar to normal distribution

🔍 OUTLIER DETECTION (IQR method):
   Lower bound: 59.207
   Upper bound: 140.760
   Potential outliers: 0 values

✅ DATA QUALITY CHECKS:
============================================================

print("EXAMPLE 1: WEBSITE TRAFFIC ANALYSIS")
print("=" * 70)

# Simulate 30 days of website traffic data
np.random.seed(123)  # For reproducible results

# Base traffic with weekly pattern (higher on weekends)
days = 30
base_traffic = 500  # Average daily visitors
weekend_boost = 150  # Extra visitors on weekends
trend = 5  # Daily growth trend

# Generate daily traffic with pattern
daily_visitors = np.zeros(days)
for day in range(days):
    # Day of week (0=Monday, 6=Sunday)
    day_of_week = day % 7

    # Base + trend + random noise
    traffic = base_traffic + (trend * day)

    # Weekend boost (Friday, Saturday, Sunday)
    if day_of_week >= 4:  # Friday=4, Saturday=5, Sunday=6
        traffic += weekend_boost

    # Random variation
    traffic += np.random.randn() * 50

    # Ensure positive
    traffic = max(traffic, 100)

    daily_visitors[day] = round(traffic)

print(f"\nGenerated {days} days of website traffic data")
print(f"First 10 days: {daily_visitors[:10].astype(int)}")

# Weekend indices (0-based, assuming Day 0 is Monday)
# In 30 days starting Monday: weekends are days 4-5, 11-12, 18-19, 25-26
weekends = np.array([4, 5, 11, 12, 18, 19, 25, 26])

print("\n" + "="*70)
print("BASIC TRAFFIC ANALYSIS")
print("="*70)

# Basic stats
total_visitors = np.sum(daily_visitors)
avg_daily = np.mean(daily_visitors)
best_day = np.max(daily_visitors)
worst_day = np.min(daily_visitors)
best_day_idx = np.argmax(daily_visitors)
worst_day_idx = np.argmin(daily_visitors)

print(f"📊 Total visitors: {total_visitors:,.0f}")
print(f"📈 Daily average: {avg_daily:.0f} visitors")
print(f"🏆 Best day: Day {best_day_idx+1} with {best_day:.0f} visitors")
print(f"📉 Worst day: Day {worst_day_idx+1} with {worst_day:.0f} visitors")
print(f"📐 Variability (std): {np.std(daily_visitors):.0f} visitors")

# Weekend vs weekday comparison
weekday_visitors = np.delete(daily_visitors, weekends)
weekend_visitors = daily_visitors[weekends]

print("\n" + "="*70)
print("WEEKDAY VS WEEKEND ANALYSIS")
print("="*70)

print(f"📅 Weekdays ({len(weekday_visitors)} days):")
print(f"   Average: {np.mean(weekday_visitors):.0f} visitors")
print(f"   Range: {np.min(weekday_visitors):.0f} - {np.max(weekday_visitors):.0f}")
print(f"   Total: {np.sum(weekday_visitors):,.0f}")

print(f"\n🎉 Weekends ({len(weekend_visitors)} days):")
print(f"   Average: {np.mean(weekend_visitors):.0f} visitors")
print(f"   Range: {np.min(weekend_visitors):.0f} - {np.max(weekend_visitors):.0f}")
print(f"   Total: {np.sum(weekend_visitors):,.0f}")

weekend_boost_pct = ((np.mean(weekend_visitors) - np.mean(weekday_visitors)) / 
                     np.mean(weekday_visitors) * 100)
print(f"\n📊 Weekend boost: +{weekend_boost_pct:.1f}% more visitors")

# Trend analysis (7-day moving average)
print("\n" + "="*70)
print("TREND ANALYSIS (7-DAY MOVING AVERAGE)")
print("="*70)

# Calculate 7-day moving average
window_size = 7
moving_avg = np.zeros(days - window_size + 1)
for i in range(len(moving_avg)):
    moving_avg[i] = np.mean(daily_visitors[i:i+window_size])

print(f"7-day moving average range: {np.min(moving_avg):.1f} - {np.max(moving_avg):.1f}")

# Check if traffic is growing
first_week_avg = np.mean(daily_visitors[:7])
last_week_avg = np.mean(daily_visitors[-7:])
growth_pct = ((last_week_avg - first_week_avg) / first_week_avg) * 100

print(f"\nFirst week average: {first_week_avg:.0f} visitors")
print(f"Last week average: {last_week_avg:.0f} visitors")
if growth_pct > 0:
    print(f"📈 Growth: +{growth_pct:.1f}% over {days} days")
else:
    print(f"📉 Decline: {growth_pct:.1f}% over {days} days")

# Performance thresholds
print("\n" + "="*70)
print("PERFORMANCE THRESHOLDS")
print("="*70)

# Define performance levels
excellent_threshold = np.percentile(daily_visitors, 75)  # Top 25%
good_threshold = np.percentile(daily_visitors, 50)       # Median
poor_threshold = np.percentile(daily_visitors, 25)       # Bottom 25%

print(f"🏆 Excellent day: > {excellent_threshold:.0f} visitors (top 25%)")
print(f"👍 Good day: {good_threshold:.0f} - {excellent_threshold:.0f} visitors")
print(f"📊 Average day: {poor_threshold:.0f} - {good_threshold:.0f} visitors")
print(f"⚠️  Poor day: < {poor_threshold:.0f} visitors (bottom 25%)")

# Count days in each category
excellent_days = np.sum(daily_visitors > excellent_threshold)
good_days = np.sum((daily_visitors <= excellent_threshold) & 
                   (daily_visitors > good_threshold))
average_days = np.sum((daily_visitors <= good_threshold) & 
                      (daily_visitors > poor_threshold))
poor_days = np.sum(daily_visitors <= poor_threshold)

print(f"\n📅 Days in each category:")
print(f"   Excellent: {excellent_days} days ({excellent_days/days*100:.0f}%)")
print(f"   Good: {good_days} days ({good_days/days*100:.0f}%)")
print(f"   Average: {average_days} days ({average_days/days*100:.0f}%)")
print(f"   Poor: {poor_days} days ({poor_days/days*100:.0f}%)")

# Outlier detection
print("\n" + "="*70)
print("OUTLIER DETECTION")
print("="*70)

q1 = np.percentile(daily_visitors, 25)
q3 = np.percentile(daily_visitors, 75)
iqr = q3 - q1
lower_bound = q1 - 1.5 * iqr
upper_bound = q3 + 1.5 * iqr

outliers = daily_visitors[(daily_visitors < lower_bound) | (daily_visitors > upper_bound)]

if len(outliers) > 0:
    print(f"🚨 Found {len(outliers)} outlier days:")
    for i, outlier in enumerate(outliers, 1):
        day_idx = np.where(daily_visitors == outlier)[0][0] + 1
        if outlier > upper_bound:
            reason = "Unusually high traffic"
        else:
            reason = "Unusually low traffic"
        print(f"   {i}. Day {day_idx}: {outlier:.0f} visitors ({reason})")
else:
    print("✅ No outlier days detected")

# Summary insights
print("\n" + "="*70)
print("KEY INSIGHTS & RECOMMENDATIONS")
print("="*70)

print("1. 📅 Weekend Strategy:") 
print(f"   - Weekends get {weekend_boost_pct:.0f}% more traffic")
print("   - Consider weekend promotions or content")

print("\n2. 📈 Growth Trend:")
if growth_pct > 5:
    print(f"   - Strong growth ({growth_pct:.0f}% increase)")
    print("   - Current strategies are working well")
elif growth_pct > 0:
    print(f"   - Moderate growth ({growth_pct:.0f}% increase)")
    print("   - Continue current efforts")
else:
    print(f"   - Traffic declined by {abs(growth_pct):.0f}%")
    print("   - Investigate causes and adjust strategy")

print("\n3. 🎯 Performance Targets:")
print(f"   - {excellent_days/days*100:.0f}% of days exceed excellent threshold")
print(f"   - Target: Increase this to 30%+")

if len(outliers) > 0:
    print("\n4. 🔍 Outlier Investigation:")
    print("   - Investigate outlier days for special events")
    print("   - Learn from high-performing days")
    print("   - Address issues on low-performing days")

print("\n" + "="*70)
print("COMPLETE STATISTICAL SUMMARY")
print("="*70)
data_summary(daily_visitors, "Website Traffic Data", round_digits=0)

print("\n\n" + "="*70)
print("EXAMPLE 2: A/B TESTING ANALYSIS")
print("="*70)

# Simulate A/B test results for a website redesign
np.random.seed(42)

# Parameters
sample_size = 1000  # Users per group
baseline_conversion = 0.12  # 12% conversion rate for old design
new_design_lift = 0.03  # 3 percentage point improvement

# Generate conversion data (1 = converted, 0 = didn't convert)
group_a_conversions = np.random.binomial(1, baseline_conversion, sample_size)
group_b_conversions = np.random.binomial(1, baseline_conversion + new_design_lift, sample_size)

print(f"\nA/B Test Setup:")
print(f"- Sample size: {sample_size:,} users per group")
print(f"- Group A (old design): Expected {baseline_conversion*100:.0f}% conversion")
print(f"- Group B (new design): Expected {(baseline_conversion + new_design_lift)*100:.0f}% conversion")
print(f"- Expected improvement: {new_design_lift*100:.0f} percentage points")

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

# Calculate conversion rates
conversions_a = np.sum(group_a_conversions)
conversions_b = np.sum(group_b_conversions)

conv_rate_a = conversions_a / sample_size * 100
conv_rate_b = conversions_b / sample_size * 100
absolute_improvement = conv_rate_b - conv_rate_a
relative_improvement = (conv_rate_b - conv_rate_a) / conv_rate_a * 100

print(f"📊 Conversion Rates:")
print(f"   Group A (old): {conv_rate_a:.2f}% ({conversions_a}/{sample_size})")
print(f"   Group B (new): {conv_rate_b:.2f}% ({conversions_b}/{sample_size})")
print(f"   Absolute improvement: {absolute_improvement:.2f} percentage points")
print(f"   Relative improvement: {relative_improvement:.1f}%")

# Statistical significance
print("\n" + "="*70)
print("STATISTICAL SIGNIFICANCE TEST")
print("="*70)

from scipy import stats

# Perform two-proportion z-test
# First, calculate pooled proportion
p1 = conversions_a / sample_size
p2 = conversions_b / sample_size
p_pool = (conversions_a + conversions_b) / (sample_size * 2)
se = np.sqrt(p_pool * (1 - p_pool) * (2 / sample_size))
z_score = (p2 - p1) / se
p_value = 2 * (1 - stats.norm.cdf(abs(z_score)))  # Two-tailed test

print(f"📈 Statistical Test Results:")
print(f"   Z-score: {z_score:.3f}")
print(f"   P-value: {p_value:.4f}")

# Common significance levels
alpha_90 = 0.10
alpha_95 = 0.05
alpha_99 = 0.01

print(f"\n🎯 Significance Levels:")
print(f"   90% confidence (α=0.10): {'Significant' if p_value < alpha_90 else 'Not significant'}")
print(f"   95% confidence (α=0.05): {'Significant' if p_value < alpha_95 else 'Not significant'}")
print(f"   99% confidence (α=0.01): {'Significant' if p_value < alpha_99 else 'Not significant'}")

# Confidence interval for the difference
print("\n" + "="*70)
print("CONFIDENCE INTERVALS")
print("="*70)

# 95% confidence interval for difference
se_diff = np.sqrt((p1 * (1 - p1) / sample_size) + (p2 * (1 - p2) / sample_size))
margin_error = 1.96 * se_diff * 100  # Convert to percentage points
ci_lower = (p2 - p1) * 100 - margin_error
ci_upper = (p2 - p1) * 100 + margin_error

print(f"95% Confidence Interval for Improvement:")
print(f"   Lower bound: {ci_lower:.2f} percentage points")
print(f"   Upper bound: {ci_upper:.2f} percentage points")
print(f"   Range: {ci_upper - ci_lower:.2f} percentage points")

# Practical significance
print("\n" + "="*70)
print("PRACTICAL SIGNIFICANCE")
print("="*70)

# Calculate business impact
monthly_visitors = 100000
average_order_value = 50
current_revenue = monthly_visitors * (conv_rate_a/100) * average_order_value
new_revenue = monthly_visitors * (conv_rate_b/100) * average_order_value
revenue_increase = new_revenue - current_revenue
annual_increase = revenue_increase * 12

print(f"💰 Business Impact (Monthly):")
print(f"   Current revenue: ${current_revenue:,.0f}")
print(f"   New design revenue: ${new_revenue:,.0f}")
print(f"   Monthly increase: ${revenue_increase:,.0f}")
print(f"   Annual increase: ${annual_increase:,.0f}")

# Minimum Detectable Effect (MDE)
print(f"\n🔍 Minimum Detectable Effect (MDE):")
# MDE at 80% power, 95% confidence
z_alpha = 1.96  # 95% confidence
z_beta = 0.84   # 80% power
mde = (z_alpha + z_beta) * np.sqrt((2 * p1 * (1 - p1)) / sample_size) * 100
print(f"   With {sample_size:,} users per group:")
print(f"   Can detect effects > {mde:.2f} percentage points")
print(f"   Our observed effect: {absolute_improvement:.2f} percentage points")
print(f"   {'✓ Effect larger than MDE' if absolute_improvement > mde else '✗ Effect smaller than MDE'}")

# Bayesian approach (simplified)
print("\n" + "="*70)
print("BAYESIAN ANALYSIS (SIMPLIFIED)")
print("="*70)

# Prior belief: 50% chance new design is better
prior_success = 0.5
prior_failure = 0.5

# Likelihood: Probability of seeing this data if new design is better
# Simplified calculation
bayes_factor = ((conv_rate_b/100) ** conversions_b * (1 - conv_rate_b/100) ** (sample_size - conversions_b)) / \
               ((conv_rate_a/100) ** conversions_a * (1 - conv_rate_a/100) ** (sample_size - conversions_a))

posterior_success = (bayes_factor * prior_success) / (bayes_factor * prior_success + prior_failure)

print(f"🤔 Bayesian Perspective:")
print(f"   Prior probability new design is better: {prior_success*100:.0f}%")
print(f"   Posterior probability: {posterior_success*100:.1f}%")
print(f"   Bayes factor: {bayes_factor:.2f} (>{10} = strong evidence)")

# Final recommendation
print("\n" + "="*70)
print("FINAL RECOMMENDATION")
print("="*70)

if p_value < alpha_95:
    if absolute_improvement > 0:
        print("✅ RECOMMENDATION: Implement the new design")
        print(f"   - Statistically significant at 95% confidence")
        print(f"   - Estimated {absolute_improvement:.1f}% improvement")
        print(f"   - Expected annual revenue increase: ${annual_increase:,.0f}")
    else:
        print("✅ RECOMMENDATION: Keep the old design")
        print(f"   - New design performed worse")
        print(f"   - {abs(absolute_improvement):.1f}% decrease in conversions")
else:
    if absolute_improvement > 0:
        print("🤔 RECOMMENDATION: Consider further testing")
        print(f"   - Not statistically significant (p={p_value:.3f})")
        print(f"   - But positive trend of {absolute_improvement:.1f}%")
        print(f"   - Run test longer or with larger sample")
    else:
        print("✅ RECOMMENDATION: Keep the old design")
        print(f"   - No evidence new design is better")
        print(f"   - {abs(absolute_improvement):.1f}% decrease in conversions")

print("\n" + "="*70)
print("DETAILED STATISTICAL SUMMARIES")
print("="*70)

print("\nGroup A (Old Design):")
data_summary(group_a_conversions, "Group A Conversions", round_digits=4)

print("\nGroup B (New Design):")
data_summary(group_b_conversions, "Group B Conversions", round_digits=4)

print("\n\n" + "="*70)
print("EXAMPLE 3: FINANCIAL DATA ANALYSIS")
print("="*70)

# Simulate stock price data
np.random.seed(101)
trading_days = 252  # Typical trading days in a year

# Parameters
initial_price = 100
daily_return_mean = 0.0005  # 0.05% daily return
daily_volatility = 0.02     # 2% daily volatility

# Generate daily returns (random walk with drift)
daily_returns = np.random.normal(daily_return_mean, daily_volatility, trading_days)

# Convert to price series
stock_prices = initial_price * np.cumprod(1 + daily_returns)

print(f"\n📊 Simulated Stock Data ({trading_days} trading days):")
print(f"- Initial price: ${initial_price:.2f}")
print(f"- Final price: ${stock_prices[-1]:.2f}")
print(f"- Daily return mean: {daily_return_mean*100:.3f}%")
print(f"- Daily volatility: {daily_volatility*100:.2f}%")

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

# Calculate returns statistics
total_return = (stock_prices[-1] / initial_price - 1) * 100
annualized_return = np.mean(daily_returns) * 252 * 100
annualized_volatility = np.std(daily_returns) * np.sqrt(252) * 100

print(f"📈 Return Metrics:")
print(f"   Total return: {total_return:.2f}%")
print(f"   Annualized return: {annualized_return:.2f}%")
print(f"   Annualized volatility: {annualized_volatility:.2f}%")

# Risk-adjusted return (Sharpe ratio)
risk_free_rate = 0.02  # 2% annual risk-free rate
daily_risk_free = risk_free_rate / 252
excess_returns = daily_returns - daily_risk_free
sharpe_ratio = np.mean(excess_returns) / np.std(excess_returns) * np.sqrt(252)

print(f"\n🎯 Risk-Adjusted Returns:")
print(f"   Sharpe ratio: {sharpe_ratio:.3f}")
print(f"   Risk-free rate assumption: {risk_free_rate*100:.1f}% annual")

# Interpret Sharpe ratio
if sharpe_ratio > 1.0:
    sharpe_eval = "Excellent"
elif sharpe_ratio > 0.5:
    sharpe_eval = "Good"
elif sharpe_ratio > 0:
    sharpe_eval = "Acceptable"
else:
    sharpe_eval = "Poor"

print(f"   Evaluation: {sharpe_eval}")

# Maximum drawdown (worst peak-to-trough decline)
print("\n" + "="*70)
print("RISK METRICS")
print("="*70)

running_max = np.maximum.accumulate(stock_prices)
drawdowns = (stock_prices - running_max) / running_max * 100
max_drawdown = np.min(drawdowns)
max_drawdown_day = np.argmin(drawdowns)

# Recovery analysis
recovery_days = None
for i in range(max_drawdown_day + 1, len(stock_prices)):
    if stock_prices[i] >= running_max[max_drawdown_day]:
        recovery_days = i - max_drawdown_day
        break

print(f"📉 Maximum Drawdown:")
print(f"   Worst decline: {abs(max_drawdown):.2f}%")
print(f"   Peak before decline: ${running_max[max_drawdown_day]:.2f}")
print(f"   Trough: ${stock_prices[max_drawdown_day]:.2f}")
if recovery_days:
    print(f"   Recovery time: {recovery_days} trading days")
else:
    print(f"   Still in drawdown at end of period")

# Value at Risk (VaR) - 95% confidence
print(f"\n💳 Value at Risk (VaR):")
var_95 = np.percentile(daily_returns, 5) * 100  # 5th percentile = worst 5% of days
print(f"   95% Daily VaR: {var_95:.2f}%")
print(f"   Interpretation: On worst 5% of days, expect loss > {abs(var_95):.2f}%")

# Expected Shortfall (CVaR)
worst_5_percent = daily_returns[daily_returns <= np.percentile(daily_returns, 5)]
expected_shortfall = np.mean(worst_5_percent) * 100
print(f"   Expected Shortfall: {expected_shortfall:.2f}%")
print(f"   Average loss on worst 5% of days: {abs(expected_shortfall):.2f}%")

# Distribution analysis
print("\n" + "="*70)
print("DISTRIBUTION ANALYSIS")
print("="*70)

from scipy import stats

# Test for normality
shapiro_stat, shapiro_p = stats.shapiro(daily_returns)
jarque_bera_stat, jarque_bera_p = stats.jarque_bera(daily_returns)

print(f"📊 Normality Tests:")
print(f"   Shapiro-Wilk test: p = {shapiro_p:.4f}")
print(f"   Jarque-Bera test: p = {jarque_bera_p:.4f}")

if shapiro_p > 0.05:
    print("   ✅ Returns appear normally distributed")
else:
    print("   ⚠️  Returns NOT normally distributed (common in finance)")

# Skewness and kurtosis
skew = stats.skew(daily_returns)
kurt = stats.kurtosis(daily_returns)

print(f"\n📈 Distribution Shape:")
print(f"   Skewness: {skew:.3f}")
if skew < -0.5:
    print("     → Left-skewed (more large negative returns)")
elif skew > 0.5:
    print("     → Right-skewed (more large positive returns)")
else:
    print("     → Approximately symmetric")

print(f"   Excess Kurtosis: {kurt:.3f}")
if kurt > 1:
    print("     → Heavy tails (more extreme returns than normal)")
elif kurt < -1:
    print("     → Light tails (fewer extreme returns than normal)")
else:
    print("     → Similar tail behavior to normal distribution")

# Autocorrelation (check for patterns)
print("\n" + "="*70)
print("SERIAL CORRELATION ANALYSIS")
print("="*70)

# Calculate autocorrelation for lags 1-5
lags = 5
autocorrelations = np.zeros(lags)
for lag in range(1, lags + 1):
    if lag < len(daily_returns):
        corr = np.corrcoef(daily_returns[:-lag], daily_returns[lag:])[0, 1]
        autocorrelations[lag-1] = corr

print(f"📅 Autocorrelation of Returns:")
for lag in range(lags):
    print(f"   Lag {lag+1}: {autocorrelations[lag]:.4f}", end="")
    if abs(autocorrelations[lag]) > 0.1:
        print(" ⚠️  (Possible pattern)")
    else:
        print(" ✅ (Random)")

# Volatility clustering test
print(f"\n📊 Volatility Clustering:")
squared_returns = daily_returns ** 2
vol_cluster_corr = np.corrcoef(squared_returns[:-1], squared_returns[1:])[0, 1]
print(f"   Correlation of squared returns: {vol_cluster_corr:.4f}")
if vol_cluster_corr > 0.1:
    print("   ⚠️  Volatility clustering present (common in financial data)")
    print("   → High volatility days tend to cluster together")
else:
    print("   ✅ No strong volatility clustering")

# Performance benchmarks
print("\n" + "="*70)
print("PERFORMANCE BENCHMARKS")
print("="*70)

# Compare to buy-and-hold
buy_hold_return = total_return

# Compare to market (simulated)
market_returns = np.random.normal(0.0003, 0.015, trading_days)  # Simpler market
market_cumulative = np.cumprod(1 + market_returns)
market_total_return = (market_cumulative[-1] / 1 - 1) * 100

# Calculate alpha and beta (simplified)
covariance = np.cov(daily_returns, market_returns)[0, 1]
market_variance = np.var(market_returns)
beta = covariance / market_variance
alpha = annualized_return - beta * (np.mean(market_returns) * 252 * 100)

print(f"📈 Performance vs Market:")
print(f"   Stock return: {total_return:.2f}%")
print(f"   Market return: {market_total_return:.2f}%")
print(f"   Outperformance: {total_return - market_total_return:.2f}%")
print(f"   Beta: {beta:.3f} (measure of market sensitivity)")
if beta > 1.2:
    print("     → High beta (more volatile than market)")
elif beta > 0.8:
    print("     → Similar volatility to market")
else:
    print("     → Low beta (less volatile than market)")

print(f"   Alpha: {alpha:.2f}% (excess return after adjusting for risk)")
if alpha > 0:
    print("     → Positive alpha (outperformed on risk-adjusted basis)")
else:
    print("     → Negative alpha (underperformed on risk-adjusted basis)")

# Monte Carlo simulation for future prices
print("\n" + "="*70)
print("MONTE CARLO SIMULATION (FORECAST)")
print("="*70)

n_simulations = 1000
forecast_days = 21  # Approximately 1 month
future_prices = np.zeros((n_simulations, forecast_days))

for i in range(n_simulations):
    # Generate future returns based on historical characteristics
    future_returns = np.random.normal(np.mean(daily_returns), 
                                      np.std(daily_returns), 
                                      forecast_days)
    future_prices[i] = stock_prices[-1] * np.cumprod(1 + future_returns)

# Analyze simulation results
final_prices = future_prices[:, -1]
median_forecast = np.median(final_prices)
percentile_5 = np.percentile(final_prices, 5)
percentile_95 = np.percentile(final_prices, 95)

print(f"🔮 1-Month Price Forecast (Monte Carlo):")
print(f"   Current price: ${stock_prices[-1]:.2f}")
print(f"   Median forecast: ${median_forecast:.2f}")
print(f"   90% Confidence interval: [${percentile_5:.2f}, ${percentile_95:.2f}]")
print(f"   Probability of gain: {(np.sum(final_prices > stock_prices[-1]) / n_simulations * 100):.1f}%")

# Risk metrics from simulation
simulation_returns = (final_prices / stock_prices[-1] - 1) * 100
var_95_sim = np.percentile(simulation_returns, 5)
expected_shortfall_sim = np.mean(simulation_returns[simulation_returns <= var_95_sim])

print(f"\n🎯 Risk Metrics from Simulation:")
print(f"   95% VaR (1-month): {var_95_sim:.2f}%")
print(f"   Expected Shortfall: {expected_shortfall_sim:.2f}%")

# Final summary
print("\n" + "="*70)
print("INVESTMENT SUMMARY & RECOMMENDATIONS")
print("="*70)

print(f"📋 Summary Statistics:")
print(f"   Annualized return: {annualized_return:.2f}%")
print(f"   Annualized volatility: {annualized_volatility:.2f}%")
print(f"   Sharpe ratio: {sharpe_ratio:.3f} ({sharpe_eval})")
print(f"   Maximum drawdown: {abs(max_drawdown):.2f}%")

print(f"\n🎯 Risk Assessment:")
if max_drawdown < -20:
    risk_level = "High"
elif max_drawdown < -10:
    risk_level = "Medium"
else:
    risk_level = "Low"

print(f"   Risk level: {risk_level}")
print(f"   Worst-case daily loss (95% VaR): {abs(var_95):.2f}%")

print(f"\n📈 Growth Potential:")
if alpha > 2:
    growth_potential = "High"
elif alpha > 0:
    growth_potential = "Moderate"
else:
    growth_potential = "Low"

print(f"   Growth potential: {growth_potential}")
print(f"   Risk-adjusted performance (alpha): {alpha:.2f}%")

print(f"\n✅ FINAL RECOMMENDATION:")
if sharpe_ratio > 0.5 and max_drawdown > -15 and alpha > 0:
    print("   CONSIDER INVESTING")
    print("   - Good risk-adjusted returns")
    print("   - Acceptable drawdown risk")
    print("   - Positive alpha indicates skill")
elif sharpe_ratio > 0:
    print("   CAUTIOUS CONSIDERATION")
    print("   - Some positive characteristics")
    print("   - Review specific risk concerns")
else:
    print("   AVOID OR REDUCE POSITION")
    print("   - Poor risk-adjusted returns")
    print("   - Consider better alternatives")

print("\n" + "="*70)
print("DETAILED STATISTICAL SUMMARIES")
print("="*70)

print("\nDaily Returns Analysis:")
data_summary(daily_returns * 100, "Daily Returns (%)", round_digits=4)

print("\nStock Price Analysis:")
data_summary(stock_prices, "Stock Prices ($)", round_digits=2)