Part 13: Simulation

import numpy as np

# create a 1-dimensional array
arr_1d = np.array([1.1, 1.4, 2.6])
print('1-Dimensional Array:', arr_1d)

# create a 2-dimensional array
arr_2d = np.array([
    [1.1, 1.4, 2.6],
    [-3.2, 5.603, 4]
])
print('\n2-Dimensional Array:')
print(arr_2d)

# slice a 1-dimensional array (like slicing a list)
print('Extract a single element:', arr_1d[1])
print('Extract a range:', arr_1d[:2])

# slice a 2-dimensional array - index is [row, column]
print('\nExtract element from 2d array:', arr_2d[0, 1])
print('Extract range from 2d array:')
print(arr_2d[:, :2])

# array addition (sum element by element)
arr_1d_second = np.array([7, 4.3, 1.1])
print('Sum of arrays:', arr_1d + arr_1d_second)

# array multiplication (multiply element by element)
print('Product of arrays:', arr_1d * arr_1d_second)

# broadcasting: multiply by scalar
print('\nBroadcast (multiply by 2):')
print(2 * arr_2d)

# sum method
print('Array:', arr_1d)
print('Sum:', arr_1d.sum())

# standard deviation method
print('\nStandard deviation:', arr_1d.std())

# other useful methods: mean(), min(), max()
print('Mean:', arr_1d.mean())

import pandas as pd
import matplotlib.pyplot as plt

# Set random seed for reproducible results
np.random.seed(42)

# Roll a six-sided die 1000 times
num_rolls = 1000
dice_rolls = np.random.randint(1, 7, num_rolls)

print(f"First 10 rolls: {dice_rolls[:10]}")
print(f"Average: {dice_rolls.mean():.2f} (expected: 3.5)")

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.hist(dice_rolls, bins=range(1, 8), alpha=0.7,
         color='skyblue', edgecolor='black')
plt.title('Distribution of Dice Rolls')
plt.xlabel('Dice Value')
plt.ylabel('Frequency')
plt.xticks(range(1, 7))

plt.subplot(1, 2, 2)
running_average = np.cumsum(dice_rolls) / np.arange(1, num_rolls + 1)
plt.plot(running_average, color='red', linewidth=2)
plt.axhline(y=3.5, color='blue', linestyle='--',
            label='Expected (3.5)')
plt.title('Running Average')
plt.xlabel('Number of Rolls')
plt.ylabel('Average Value')
plt.legend()
plt.tight_layout()

# Monte Carlo simulation: monthly sales forecast

base_monthly_sales = 25000
monthly_variation = 0.20  # 20% variation

# Seasonal factors for each month
seasonal_factors = [
    0.9, 0.9, 1.0, 1.0, 1.1, 1.1,  # Jan-Jun
    1.0, 1.0, 1.0, 1.1, 1.2, 1.3   # Jul-Dec
]

def simulate_one_year():
    monthly_sales = []
    for month in range(12):
        expected_sales = base_monthly_sales * seasonal_factors[month]
        actual_sales = np.random.normal(expected_sales,
                                       expected_sales * monthly_variation)
        actual_sales = max(actual_sales, 0)
        monthly_sales.append(actual_sales)
    return monthly_sales

# Simulate 1,000 possible years
num_simulations = 1000
all_annual_totals = []

for sim in range(num_simulations):
    year_sales = simulate_one_year()
    all_annual_totals.append(sum(year_sales))

all_annual_totals = np.array(all_annual_totals)

print(f"Average annual sales: ${all_annual_totals.mean():,.0f}")
print(f"Minimum: ${all_annual_totals.min():,.0f}")
print(f"Maximum: ${all_annual_totals.max():,.0f}")
print(f"10th percentile: ${np.percentile(all_annual_totals, 10):,.0f}")
print(f"90th percentile: ${np.percentile(all_annual_totals, 90):,.0f}")

plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.hist(all_annual_totals, bins=50, alpha=0.7,
         color='lightblue', edgecolor='black')
plt.axvline(all_annual_totals.mean(), color='red',
            linestyle='--', label=f'Mean: ${all_annual_totals.mean():,.0f}')
plt.title('Distribution of Annual Sales')
plt.xlabel('Annual Sales ($)')
plt.ylabel('Frequency')
plt.legend()

plt.subplot(1, 2, 2)
month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
               'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
example_year = simulate_one_year()
plt.bar(month_names, example_year, color='lightgreen', alpha=0.7)
plt.title('One Example Year')
plt.xlabel('Month')
plt.ylabel('Sales ($)')
plt.xticks(rotation=45)
plt.tight_layout()

# Business question: What's the probability of meeting target?
target_sales = 300000
prob_meeting_target = np.mean(all_annual_totals >= target_sales) * 100

print(f"Probability of meeting ${target_sales:,} target: {prob_meeting_target:.1f}%")

# Show different targets
for target in [280000, 300000, 320000]:
    prob = np.mean(all_annual_totals >= target) * 100
    print(f"Target ${target:,}: {prob:.1f}% probability")