import numpy as np
from scipy.stats import binom, beta

np.random.seed(42)
N = 5000  # simulation replications


# Accident probability & trips
p_acc = 1/1000
n_trips = 500

# Beta cost distribution parameters
alpha, beta_param = 1.5, 3
cost_min, cost_max = 500, 20000

total_costs = []

for _ in range(N):
    # Number of accidents in 500 trips
    num_acc = binom.rvs(n_trips, p_acc)

    # Cost of each accident if any occur
    if num_acc > 0:
        raw_costs = beta.rvs(alpha, beta_param, size=num_acc)
        scaled_costs = cost_min + raw_costs * (cost_max - cost_min)
        total_costs.append(np.sum(scaled_costs))
    else:
        total_costs.append(0)

total_costs = np.array(total_costs)

# Mean + CI
mean_cost = total_costs.mean()
std_cost = total_costs.std()
ci_mean = (
    mean_cost - 1.96 * std_cost / np.sqrt(N),
    mean_cost + 1.96 * std_cost / np.sqrt(N)
)

# Probability total cost > 8000
indicator = (total_costs > 8000).astype(int)
p_exceed = indicator.mean()
ci_prop = (
    p_exceed - 1.96*np.sqrt(p_exceed*(1-p_exceed)/N),
    p_exceed + 1.96*np.sqrt(p_exceed*(1-p_exceed)/N)
)

# PRINT RESULTS
print("\n1.1 Type of random variable for # of accidents:")
print("    → Binomial(n = 500, p = 0.001)\n")

print("1.3 Average total accident cost for 500 trips:")
print(f"    Mean total cost = ${mean_cost:,.2f}")
print(f"    95% CI for mean cost = (${ci_mean[0]:,.2f}, ${ci_mean[1]:,.2f})\n")

print("1.4 Probability that total accident cost exceeds $8,000:")
print(f"    Estimated probability = {p_exceed:.4f}")
print(f"    95% CI for probability = ({ci_prop[0]:.4f}, {ci_prop[1]:.4f})\n")

print("Interpretation:")
print("    The mean cost represents the expected financial risk over 500 trips.")
print("    The confidence interval expresses uncertainty about this estimate.")
print("    The exceedance probability tells how often a high-cost scenario (> $8,000) occurs.\n")

1.1 Type of random variable for # of accidents:
    → Binomial(n = 500, p = 0.001)

1.3 Average total accident cost for 500 trips:
    Mean total cost = $3,592.14
    95% CI for mean cost = ($3,432.54, $3,751.74)

1.4 Probability that total accident cost exceeds $8,000:
    Estimated probability = 0.1998
    95% CI for probability = (0.1887, 0.2109)

Interpretation:
    The mean cost represents the expected financial risk over 500 trips.
    The confidence interval expresses uncertainty about this estimate.
    The exceedance probability tells how often a high-cost scenario (> $8,000) occurs.

import numpy as np

np.random.seed(42)
N = 5000



# Activity Distributions
A_times=[5,6,7,8]; A_probs=[0.25,0.35,0.25,0.15]
B_times=[3,5,7]; B_probs=[0.20,0.55,0.25]
C_times=[10,12,14,16,18]; C_probs=[0.10,0.25,0.40,0.20,0.05]
D_times=[8,10]; D_probs=[0.60,0.40]

project_lengths=[]

for _ in range(N):
    A = np.random.choice(A_times, p=A_probs)
    B = np.random.choice(B_times, p=B_probs)
    C = np.random.choice(C_times, p=C_probs)
    D = np.random.choice(D_times, p=D_probs)
    project_lengths.append(A + B + C + D)

project_lengths = np.array(project_lengths)

# Mean & Std Dev
mean_T = project_lengths.mean()
std_T = project_lengths.std()

# 95% CI for Mean
ci_mean = (
    mean_T - 1.96 * std_T / np.sqrt(N),
    mean_T + 1.96 * std_T / np.sqrt(N)
)

# Probability of finishing ≤ 35 weeks
prob_35 = np.mean(project_lengths <= 35)

# 95% CI for Probability
ci_prob = (
    prob_35 - 1.96 * np.sqrt(prob_35*(1-prob_35)/N),
    prob_35 + 1.96 * np.sqrt(prob_35*(1-prob_35)/N)
)

# PRINT RESULTS
print("\n2.1 Project Duration Statistics:")
print(f"    Average project length (weeks) = {mean_T:.4f}")
print(f"    Standard deviation = {std_T:.4f}")
print(f"    95% CI for mean duration = ({ci_mean[0]:.4f}, {ci_mean[1]:.4f})\n")

print("2.2 Probability of finishing in 35 weeks or less:")
print(f"    Estimated probability = {prob_35:.4f}")
print(f"    95% CI for this probability = ({ci_prob[0]:.4f}, {ci_prob[1]:.4f})\n")

print("Interpretation:")
print("    The 95% confidence interval for the mean project duration shows the range")
print("    in which the true average duration is likely to fall.")
print("    The probability confidence interval quantifies uncertainty in meeting the")
print("    35-week deadline, helping managers evaluate schedule risk.\n")

2.1 Project Duration Statistics:
    Average project length (weeks) = 33.8904
    Standard deviation = 2.8221
    95% CI for mean duration = (33.8122, 33.9686)

2.2 Probability of finishing in 35 weeks or less:
    Estimated probability = 0.7144
    95% CI for this probability = (0.7019, 0.7269)

Interpretation:
    The 95% confidence interval for the mean project duration shows the range
    in which the true average duration is likely to fall.
    The probability confidence interval quantifies uncertainty in meeting the
    35-week deadline, helping managers evaluate schedule risk.

import numpy as np

np.random.seed(42)
N = 5000

price = 200
fixed_cost = 350000

profits = []

for _ in range(N):
    vc = np.random.uniform(90,150)         # variable cost per unit
    D  = max(0, np.random.normal(9000,2000))  # demand truncated at zero
    profit = (price - vc)*D - fixed_cost
    profits.append(profit)

profits = np.array(profits)

# Mean profit
mean_profit = profits.mean()
std_profit  = profits.std()

ci_mean = (
    mean_profit - 1.96 * std_profit / np.sqrt(N),
    mean_profit + 1.96 * std_profit / np.sqrt(N)
)

# Probability of loss
p_loss = np.mean(profits < 0)
ci_loss = (
    p_loss - 1.96 * np.sqrt(p_loss*(1-p_loss)/N),
    p_loss + 1.96 * np.sqrt(p_loss*(1-p_loss)/N)
)

print("\nProblem 3 Results:")
print(f"Average Profit = ${mean_profit:,.2f}")
print(f"95% CI for Mean Profit = (${ci_mean[0]:,.2f}, ${ci_mean[1]:,.2f})")

print(f"\nProbability of Loss = {p_loss:.4f}")
print(f"95% CI for Probability of Loss = ({ci_loss[0]:.4f}, {ci_loss[1]:.4f})")

Problem 3:

Average Profit = $372,838.70
95% CI for Mean Profit = ($366,554.15, $379,123.24)

Probability of Loss = 0.0260
95% CI for Probability of Loss = (0.0216, 0.0304)