import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import pulp as pl
from itertools import combinations
import matplotlib.cm as cm
import matplotlib.colors as mcolors

# Set random seed for reproducibility
np.random.seed(42)

# Parameters
num_vertices = 20
edge_probability = 0.2  # Probability of an edge between any two vertices

# Generate a random graph
G = nx.gnp_random_graph(num_vertices, edge_probability)

# We can also add some additional edges to make it more challenging
# Add edges between some specific nodes to create more constraints
additional_edges = [(0, 5), (1, 6), (2, 7), (3, 8), (4, 9)]
for u, v in additional_edges:
    if u < num_vertices and v < num_vertices and not G.has_edge(u, v):
        G.add_edge(u, v)

# Print basic graph statistics
print(f"Graph with {G.number_of_nodes()} vertices and {G.number_of_edges()} edges")
print(f"Average degree: {2 * G.number_of_edges() / G.number_of_nodes():.2f}")

# Upper bound on the number of colors needed (worst case)
max_colors = min(num_vertices, max(dict(G.degree()).values()) + 1)
print(f"Maximum number of colors needed (worst case): {max_colors}")

# Create the optimization model
model = pl.LpProblem("Graph_Coloring", pl.LpMinimize)

# Set of colors (using indices)
colors = range(max_colors)

# Decision variables
# x[i, v] = 1 if vertex i is assigned color v, 0 otherwise
x = {(i, v): pl.LpVariable(f"x_{i}_{v}", cat='Binary') 
     for i in G.nodes() for v in colors}

# y[v] = 1 if color v is used, 0 otherwise
y = {v: pl.LpVariable(f"y_{v}", cat='Binary') for v in colors}

# Objective function: minimize the number of colors used
model += pl.lpSum(y[v] for v in colors)

# Constraint: each vertex must be assigned exactly one color
for i in G.nodes():
    model += pl.lpSum(x[(i, v)] for v in colors) == 1

# Constraint: adjacent vertices must have different colors
for (i, j) in G.edges():
    for v in colors:
        model += x[(i, v)] + x[(j, v)] <= 1

# Constraint: if a vertex is assigned a color, that color is used
for i in G.nodes():
    for v in colors:
        model += x[(i, v)] <= y[v]

# Constraint: ensure colors are used in order (symmetry breaking)
for v in range(1, max_colors):
    model += y[v] <= y[v-1]

# Solve the model
model.solve(pl.PULP_CBC_CMD(msg=False))

# Extract results
vertex_colors = {}
colors_used = set()

for i in G.nodes():
    for v in colors:
        if pl.value(x[(i, v)]) > 0.5:
            vertex_colors[i] = v
            colors_used.add(v)

num_colors_used = len(colors_used)

print(f"Optimal solution status: {pl.LpStatus[model.status]}")
print(f"Minimum number of colors needed: {num_colors_used}")
print("Vertex coloring:")
for i, color in sorted(vertex_colors.items()):
    print(f"Vertex {i}: Color {color}")

# Compare with greedy coloring
greedy_coloring = nx.greedy_color(G)
num_greedy_colors = max(greedy_coloring.values()) + 1
print(f"Number of colors with greedy algorithm: {num_greedy_colors}")

# Heuristic: DSatur algorithm
dsatur_coloring = nx.coloring.greedy_color(G, strategy="saturation_largest_first")
num_dsatur_colors = max(dsatur_coloring.values()) + 1
print(f"Number of colors with DSatur algorithm: {num_dsatur_colors}")

# Visualization of the colored graph
plt.figure(figsize=(12, 10))

# Choose a color map for visualization
color_names = list(mcolors.TABLEAU_COLORS.keys())
selected_colors = color_names[:max(num_colors_used, num_greedy_colors, num_dsatur_colors)]

# Create position layout for nodes
pos = nx.spring_layout(G, seed=42)

# Create a subplot for the optimal solution
plt.subplot(2, 2, 1)
node_colors_optimal = [vertex_colors[i] for i in G.nodes()]
color_map_optimal = [selected_colors[c] for c in node_colors_optimal]

nx.draw(G, pos, with_labels=True, node_color=color_map_optimal, 
        node_size=500, font_color="white", font_weight="bold",
        edge_color="gray", width=2, alpha=0.9)

plt.title(f"Optimal Coloring: {num_colors_used} Colors")

# Create a subplot for the greedy solution
plt.subplot(2, 2, 2)
node_colors_greedy = [greedy_coloring[i] for i in G.nodes()]
color_map_greedy = [selected_colors[c] for c in node_colors_greedy]

nx.draw(G, pos, with_labels=True, node_color=color_map_greedy, 
        node_size=500, font_color="white", font_weight="bold",
        edge_color="gray", width=2, alpha=0.9)

plt.title(f"Greedy Coloring: {num_greedy_colors} Colors")

# Create a subplot for the DSatur solution
plt.subplot(2, 2, 3)
node_colors_dsatur = [dsatur_coloring[i] for i in G.nodes()]
color_map_dsatur = [selected_colors[c] for c in node_colors_dsatur]

nx.draw(G, pos, with_labels=True, node_color=color_map_dsatur, 
        node_size=500, font_color="white", font_weight="bold",
        edge_color="gray", width=2, alpha=0.9)

plt.title(f"DSatur Coloring: {num_dsatur_colors} Colors")

# Generate a legend for the colors
plt.subplot(2, 2, 4)
plt.axis('off')
for i in range(max(num_colors_used, num_greedy_colors, num_dsatur_colors)):
    plt.plot([0], [0], 'o', color=selected_colors[i], label=f'Color {i}')
plt.legend(loc='center', fontsize=12)
plt.title("Color Legend")

plt.tight_layout()

# Algorithm performance comparison
plt.figure(figsize=(10, 6))
methods = ['Optimal', 'Greedy', 'DSatur']
colors_count = [num_colors_used, num_greedy_colors, num_dsatur_colors]

bars = plt.bar(methods, colors_count, color=['green', 'blue', 'orange'])
plt.ylabel('Number of Colors Used')
plt.title('Comparison of Graph Coloring Algorithms')
plt.grid(axis='y', linestyle='--', alpha=0.7)

# Add value labels on top of bars
for bar in bars:
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height + 0.05,
            f'{height}', ha='center', va='bottom')

# Graph statistics visualization
plt.figure(figsize=(12, 6))

# Calculate degree distribution
degree_sequence = sorted([d for n, d in G.degree()], reverse=True)
degree_count = {}
for d in degree_sequence:
    if d in degree_count:
        degree_count[d] += 1
    else:
        degree_count[d] = 1

plt.subplot(1, 2, 1)
plt.bar(degree_count.keys(), degree_count.values(), color='skyblue')
plt.title('Degree Distribution')
plt.xlabel('Degree')
plt.ylabel('Count')
plt.grid(axis='y', linestyle='--', alpha=0.7)

# Create conflict matrix visualization
plt.subplot(1, 2, 2)
adj_matrix = nx.to_numpy_array(G)
plt.imshow(adj_matrix, cmap='Blues', interpolation='none')
plt.colorbar(label='Connection')
plt.title('Adjacency Matrix')
plt.xlabel('Vertex')
plt.ylabel('Vertex')

plt.tight_layout()

# Chromatic number analysis based on graph density
plt.figure(figsize=(10, 6))

# Calculate graph density
density = nx.density(G)
print(f"Graph density: {density:.3f}")

# Theoretical bounds
vertices = G.number_of_nodes()
edges = G.number_of_edges()
max_degree = max(dict(G.degree()).values())
brooks_bound = max(2, min(vertices, max_degree))

# Generate points for the theoretical bounds
densities = np.linspace(0, 1, 100)
clique_bounds = [max(2, int(vertices * d) + 1) for d in densities]  # Rough estimate based on clique number

plt.plot(densities, clique_bounds, 'r--', label='Theoretical upper bound')
plt.scatter([density], [num_colors_used], color='green', s=100, label='Optimal solution')
plt.scatter([density], [num_greedy_colors], color='blue', s=100, label='Greedy algorithm')
plt.scatter([density], [num_dsatur_colors], color='orange', s=100, label='DSatur algorithm')
plt.scatter([density], [brooks_bound], color='purple', s=100, label="Brooks' bound")

plt.axvline(x=density, color='gray', linestyle='--', alpha=0.5)
plt.text(density+0.02, 2, f'Current density: {density:.3f}', rotation=90, alpha=0.7)

plt.xlabel('Graph Density')
plt.ylabel('Chromatic Number (Colors Required)')
plt.title('Relationship Between Graph Density and Chromatic Number')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.linalg import expm
from scipy.special import hermite, factorial

# Set plotting parameters for better visualization
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['font.size'] = 14

# Physical constants
hbar = 1.0  # Reduced Planck constant (using natural units)
m = 1.0     # Mass of the particle
omega = 1.0 # Angular frequency of the oscillator

# Simulation parameters
n_points = 1000  # Number of position points for discretization
x_range = 10.0   # Position range for simulation
dx = 2.0 * x_range / n_points  # Position step size
x = np.linspace(-x_range, x_range, n_points)  # Position array

# Create potential energy function (harmonic oscillator: V(x) = 0.5 * m * omega^2 * x^2)
V = 0.5 * m * omega**2 * x**2

# Create second derivative operator (kinetic energy part)
diag = np.ones(n_points) * (2.0 / dx**2)
off_diag = np.ones(n_points-1) * (-1.0 / dx**2)
D2 = sparse.diags([off_diag, diag, off_diag], [-1, 0, 1])

# Create Hamiltonian operator: H = -hbar^2/(2m) * d^2/dx^2 + V(x)
H = sparse.diags(V) - (hbar**2 / (2.0 * m)) * D2

# Find eigenvalues and eigenvectors of the Hamiltonian
eigenvalues, eigenvectors = sparse.linalg.eigs(H, k=10, which='SM')  # Find lowest 10 energy states
idx = eigenvalues.argsort()  # Sort by energy
eigenvalues = eigenvalues[idx].real
eigenvectors = eigenvectors[:, idx]

# Normalize eigenvectors
for i in range(len(eigenvalues)):
    eigenvectors[:, i] = eigenvectors[:, i] / np.sqrt(np.sum(np.abs(eigenvectors[:, i])**2) * dx)

# Analytical solutions for comparison
def analytical_wavefunction(n, x):
    """
    Calculate the analytical wavefunction for quantum harmonic oscillator
    for quantum number n at position x
    """
    # Constants for the harmonic oscillator
    alpha = m * omega / hbar
    
    # Normalization constant
    normalization = (alpha / np.pi)**0.25 * (1.0 / np.sqrt(2**n * factorial(n)))
    
    # Hermite polynomial
    herm = hermite(n)(np.sqrt(alpha) * x)
    
    # Full wavefunction
    psi = normalization * herm * np.exp(-alpha * x**2 / 2.0)
    
    return psi

# Calculate analytical energy levels
def analytical_energy(n):
    """Calculate analytical energy for quantum harmonic oscillator with quantum number n"""
    return hbar * omega * (n + 0.5)

# Compute analytical solutions
analytical_energies = [analytical_energy(n) for n in range(10)]
analytical_wavefunctions = [analytical_wavefunction(n, x) for n in range(10)]

# Plot and compare numerical and analytical solutions
fig, axs = plt.subplots(5, 2, figsize=(15, 20))
axs = axs.flatten()

for i in range(10):
    axs[i].plot(x, eigenvectors[:, i].real, 'b-', label=f'Numerical ψ{i}')
    axs[i].plot(x, analytical_wavefunctions[i], 'r--', label=f'Analytical ψ{i}')
    axs[i].axhline(y=0, color='k', linestyle='-', alpha=0.3)
    axs[i].set_title(f'Energy Level {i}: E_numerical = {eigenvalues[i]:.5f}, E_analytical = {analytical_energies[i]:.5f}')
    axs[i].set_xlabel('Position (x)')
    axs[i].set_ylabel('Wave Function (ψ)')
    axs[i].legend()
    axs[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Time evolution of a wavepacket
# Create a Gaussian wavepacket as initial state (superposition of ground and first excited states)
psi_0 = 0.7 * eigenvectors[:, 0] + 0.7 * eigenvectors[:, 1]
psi_0 = psi_0 / np.sqrt(np.sum(np.abs(psi_0)**2) * dx)  # Normalize

# Time parameters
t_total = 10.0  # Total simulation time
n_steps = 100   # Number of time steps
dt = t_total / n_steps  # Time step

# Store wavefunction at different times
psi_t = np.zeros((n_steps, n_points), dtype=complex)
psi_t[0] = psi_0

# Time evolution operator
def time_evolution(psi, H, dt):
    """Time evolution using the operator exp(-i*H*dt/hbar)"""
    # We use the fact that eigenvectors of H provide a basis for expansion
    coeffs = np.zeros(10, dtype=complex)
    for i in range(10):
        coeffs[i] = np.sum(np.conj(eigenvectors[:, i]) * psi) * dx
    
    # Evolve each component and recombine
    psi_evolved = np.zeros_like(psi, dtype=complex)
    for i in range(10):
        psi_evolved += coeffs[i] * np.exp(-1j * eigenvalues[i] * dt / hbar) * eigenvectors[:, i]
    
    return psi_evolved

# Perform time evolution
for t in range(1, n_steps):
    psi_t[t] = time_evolution(psi_t[t-1], H, dt)

# Probability densities for animation
prob_densities = np.abs(psi_t)**2

# Plot time evolution
fig, ax = plt.subplots(figsize=(12, 8))

# Plot potential energy curve (scaled for visibility)
V_scaled = V / np.max(V) * np.max(prob_densities) * 0.8
ax.plot(x, V_scaled, 'k-', alpha=0.3, label='Potential (scaled)')

# Select specific timesteps to plot
time_points = [0, 25, 50, 75, 99]
colors = ['b', 'g', 'r', 'c', 'm']

for i, t in enumerate(time_points):
    ax.plot(x, prob_densities[t], colors[i], label=f't = {t*dt:.2f}')

ax.set_xlabel('Position (x)')
ax.set_ylabel('Probability Density |ψ(x,t)|²')
ax.set_title('Time Evolution of Quantum Harmonic Oscillator Wavepacket')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Create an energy level diagram
n_levels = 10
fig, ax = plt.subplots(figsize=(8, 10))

# Plot energy levels
for n in range(n_levels):
    E = analytical_energy(n)
    ax.plot([-0.5, 0.5], [E, E], 'b-', linewidth=2)
    ax.text(0.6, E, f'n = {n}, E = {E:.2f}ħω', va='center')

# Add potential well visualization
x_pot = np.linspace(-0.4, 0.4, 100)
V_pot = 0.5 * m * omega**2 * x_pot**2
ax.plot(x_pot, V_pot, 'r--', label='V(x) = ½mω²x²')

ax.set_xlim(-1, 1.5)
ax.set_ylim(-0.5, analytical_energy(n_levels-1) + 1)
ax.set_ylabel('Energy (in units of ħω)')
ax.set_title('Energy Levels of Quantum Harmonic Oscillator')
ax.set_xticks([])
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

# Set random seed for reproducibility
np.random.seed(42)

# Parameters for our stock price model
initial_price = 100
drift = 0.0002  # Daily drift (expected return)
volatility = 0.015  # Daily volatility
trading_days = 252  # Number of trading days in a year
num_simulations = 1000  # Number of simulations to run
time_horizon = 252  # Simulate for 1 year

# Create an array to store all simulation results
all_simulations = np.zeros((num_simulations, time_horizon + 1))
all_simulations[:, 0] = initial_price  # Set initial price

# Perform the simulations using Geometric Brownian Motion
for i in range(num_simulations):
    # Generate random shocks
    random_shocks = np.random.normal(0, 1, time_horizon)
    
    # Calculate daily returns
    daily_returns = np.exp(drift - 0.5 * volatility**2 + volatility * random_shocks)
    
    # Calculate price path
    price_path = initial_price * np.cumprod(daily_returns)
    
    # Store the simulation
    all_simulations[i, 1:] = price_path

# Calculate statistics
final_prices = all_simulations[:, -1]
mean_final_price = np.mean(final_prices)
median_final_price = np.median(final_prices)
std_final_price = np.std(final_prices)

# Calculate probability of ending above certain thresholds
prob_above_initial = np.mean(final_prices > initial_price)
prob_above_110 = np.mean(final_prices > 110)
prob_above_120 = np.mean(final_prices > 120)
prob_below_90 = np.mean(final_prices < 90)

# Calculate confidence intervals
confidence_level = 0.95
lower_percentile = (1 - confidence_level) / 2
upper_percentile = 1 - lower_percentile
confidence_interval = np.percentile(final_prices, [lower_percentile * 100, upper_percentile * 100])

# Plotting
plt.figure(figsize=(15, 10))

# Plot 1: Sample paths
plt.subplot(2, 2, 1)
for i in range(50):  # Plot just 50 sample paths for clarity
    plt.plot(all_simulations[i], linewidth=0.8, alpha=0.6)
plt.title('Sample Stock Price Paths')
plt.xlabel('Trading Days')
plt.ylabel('Stock Price')
plt.grid(True)

# Plot 2: Histogram of final prices
plt.subplot(2, 2, 2)
sns.histplot(final_prices, kde=True, bins=30)
plt.axvline(initial_price, color='r', linestyle='--', label='Initial Price')
plt.axvline(mean_final_price, color='g', linestyle='-', label='Mean Final Price')
plt.axvline(confidence_interval[0], color='orange', linestyle=':', label=f'{confidence_level*100}% CI Lower')
plt.axvline(confidence_interval[1], color='orange', linestyle=':', label=f'{confidence_level*100}% CI Upper')
plt.title('Distribution of Final Stock Prices')
plt.xlabel('Stock Price')
plt.ylabel('Frequency')
plt.legend()
plt.grid(True)

# Plot 3: Empirical CDF
plt.subplot(2, 2, 3)
x = np.sort(final_prices)
y = np.arange(1, len(x) + 1) / len(x)
plt.plot(x, y)
plt.axvline(initial_price, color='r', linestyle='--', label='Initial Price')
plt.axhline(prob_above_initial, color='b', linestyle='--')
plt.axvline(110, color='g', linestyle='--', label='$110')
plt.axhline(prob_above_110, color='b', linestyle=':')
plt.title('Empirical Cumulative Distribution Function')
plt.xlabel('Stock Price')
plt.ylabel('Probability')
plt.grid(True)
plt.legend()

# Plot 4: Q-Q plot to check for normality of returns
plt.subplot(2, 2, 4)
log_returns = np.log(all_simulations[:, 1:] / all_simulations[:, :-1])
log_returns = log_returns.flatten()
stats.probplot(log_returns, dist="norm", plot=plt)
plt.title('Q-Q Plot of Log Returns')
plt.grid(True)

plt.tight_layout()
plt.savefig('stock_price_simulation.png', dpi=300)
plt.show()

# Display the statistical results
print(f"Initial stock price: ${initial_price:.2f}")
print(f"Mean final price: ${mean_final_price:.2f}")
print(f"Median final price: ${median_final_price:.2f}")
print(f"Standard deviation of final price: ${std_final_price:.2f}")
print(f"95% confidence interval: [${confidence_interval[0]:.2f}, ${confidence_interval[1]:.2f}]")
print(f"Probability of ending above initial price (${initial_price:.2f}): {prob_above_initial:.2%}")
print(f"Probability of ending above $110: {prob_above_110:.2%}")
print(f"Probability of ending above $120: {prob_above_120:.2%}")
print(f"Probability of ending below $90: {prob_below_90:.2%}")

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from matplotlib import cm
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D

# Black-Scholes Option Pricing Model
def black_scholes(S, K, T, r, sigma, option='call'):
    """
    Calculate Black-Scholes option price for a call or put
    
    Parameters:
    S: Current stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        price = S * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        price = K * np.exp(-r * T) * stats.norm.cdf(-d2) - S * stats.norm.cdf(-d1)
        
    return price

# Example 1: Visualize call option price as a function of stock price and time to maturity
# Parameters
S_range = np.linspace(80, 120, 50)  # Stock price range
T_range = np.linspace(0.1, 1, 50)   # Time to maturity range
K = 100                             # Strike price
r = 0.05                            # Risk-free rate
sigma = 0.2                         # Volatility

# Create meshgrid
S_mesh, T_mesh = np.meshgrid(S_range, T_range)
option_price = np.zeros_like(S_mesh)

# Calculate option price for each combination
for i in range(len(T_range)):
    for j in range(len(S_range)):
        option_price[i, j] = black_scholes(S_mesh[i, j], K, T_mesh[i, j], r, sigma, 'call')

# Plot the 3D surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(S_mesh, T_mesh, option_price, cmap=cm.coolwarm, alpha=0.8)

# Add labels
ax.set_xlabel('Stock Price ($)')
ax.set_ylabel('Time to Maturity (years)')
ax.set_zlabel('Call Option Price ($)')
ax.set_title('Black-Scholes Call Option Price')

# Add a color bar
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

plt.savefig('black_scholes_3d.png')
plt.show()

# Example 2: Visualize the Greeks (Delta, Gamma, Theta)
def bs_delta(S, K, T, r, sigma, option='call'):
    """Calculate option Delta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        return stats.norm.cdf(d1)
    else:
        return stats.norm.cdf(d1) - 1

def bs_gamma(S, K, T, r, sigma):
    """Calculate option Gamma (same for calls and puts)"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    return stats.norm.pdf(d1) / (S * sigma * np.sqrt(T))

def bs_theta(S, K, T, r, sigma, option='call'):
    """Calculate option Theta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * stats.norm.cdf(-d2)
        
    return theta / 365  # Convert to daily

# Calculate Greeks for various stock prices
S_range = np.linspace(80, 120, 100)
T = 0.5  # 6 months

delta_call = [bs_delta(S, K, T, r, sigma, 'call') for S in S_range]
delta_put = [bs_delta(S, K, T, r, sigma, 'put') for S in S_range]
gamma = [bs_gamma(S, K, T, r, sigma) for S in S_range]
theta_call = [bs_theta(S, K, T, r, sigma, 'call') for S in S_range]
theta_put = [bs_theta(S, K, T, r, sigma, 'put') for S in S_range]

# Plot the Greeks
fig, axs = plt.subplots(3, 1, figsize=(12, 15))

# Delta
axs[0].plot(S_range, delta_call, 'b-', label='Call Delta')
axs[0].plot(S_range, delta_put, 'r-', label='Put Delta')
axs[0].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[0].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[0].set_xlabel('Stock Price ($)')
axs[0].set_ylabel('Delta')
axs[0].set_title('Option Delta vs Stock Price')
axs[0].grid(True)
axs[0].legend()

# Gamma
axs[1].plot(S_range, gamma, 'g-')
axs[1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[1].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[1].set_xlabel('Stock Price ($)')
axs[1].set_ylabel('Gamma')
axs[1].set_title('Option Gamma vs Stock Price')
axs[1].grid(True)

# Theta
axs[2].plot(S_range, theta_call, 'b-', label='Call Theta')
axs[2].plot(S_range, theta_put, 'r-', label='Put Theta')
axs[2].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[2].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[2].set_xlabel('Stock Price ($)')
axs[2].set_ylabel('Theta ($ per day)')
axs[2].set_title('Option Theta vs Stock Price')
axs[2].grid(True)
axs[2].legend()

plt.tight_layout()
plt.savefig('option_greeks.png')
plt.show()

# Example 3: Monte Carlo simulation for option pricing
def monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations=10000, option='call'):
    """
    Price options using Monte Carlo simulation
    
    Parameters:
    S0: Initial stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    num_simulations: Number of simulation paths
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    # Generate random simulations of stock price at expiration
    z = np.random.standard_normal(num_simulations)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)
    
    # Calculate payoff at expiration
    if option == 'call':
        payoff = np.maximum(ST - K, 0)
    else:
        payoff = np.maximum(K - ST, 0)
    
    # Calculate present value
    option_price = np.exp(-r * T) * np.mean(payoff)
    
    return option_price, ST

# Run Monte Carlo simulation
np.random.seed(42)  # For reproducibility
S0 = 100
num_simulations = 10000

mc_price, ST = monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations, 'call')
bs_price = black_scholes(S0, K, T, r, sigma, 'call')

# Create histogram of simulated stock prices at expiration
plt.figure(figsize=(12, 6))
plt.hist(ST, bins=50, alpha=0.5, color='blue', density=True)

# Add a normal distribution curve for comparison
x = np.linspace(min(ST), max(ST), 100)
mu = S0 * np.exp(r * T)  # Expected value under risk-neutral measure
std = S0 * np.exp(r * T) * np.sqrt(np.exp(sigma**2 * T) - 1)
plt.plot(x, stats.norm.pdf(x, mu, std), 'r-', linewidth=2)

plt.axvline(x=K, color='k', linestyle='--', label=f'Strike Price (K={K})')
plt.axvline(x=S0, color='g', linestyle='-', label=f'Initial Price (S0={S0})')

plt.title(f'Monte Carlo Simulation of Stock Prices at Expiration (T={T} years)')
plt.xlabel('Stock Price at Expiration ($)')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.savefig('monte_carlo_simulation.png')
plt.show()

# Compare Monte Carlo with Black-Scholes
print(f"Monte Carlo Call Option Price: ${mc_price:.4f}")
print(f"Black-Scholes Call Option Price: ${bs_price:.4f}")
print(f"Difference: ${abs(mc_price - bs_price):.4f}")

# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
from sympy import primepi
from math import log

# Define the range and step size
x_values = np.linspace(10, 100000, 200)  # Range from 10 to 100,000 with 200 points

# Calculate prime counting function pi(x)
pi_values = [primepi(int(x)) for x in x_values]

# Calculate the Prime Number Theorem (PNT) approximation
pnt_approx = [x / log(x) for x in x_values]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x_values, pi_values, label=r'$\pi(x)$ (Prime Count)', color='blue')
plt.plot(x_values, pnt_approx, label=r'$\frac{x}{\log x}$ (PNT Approximation)', linestyle='--', color='red')

# Add title and labels
plt.title("Prime Number Theorem: $\pi(x)$ vs. $\dfrac{x}{\log x}$", fontsize=14)
plt.xlabel("x", fontsize=12)
plt.ylabel("Value", fontsize=12)
plt.legend()
plt.grid(True)

# Show the plot
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d import Axes3D

# Define the linear transformation matrix
T = np.array([[2, 1], 
              [1, 3]])

print("Linear Transformation Matrix T:")
print(T)

# Part 1: Transform some example points
points = np.array([[1, 0], [0, 1], [1, 1], [-1, 2]])
transformed_points = np.dot(points, T)

print("\nOriginal Points:")
for p in points:
    print(f"[{p[0]}, {p[1]}]")
    
print("\nTransformed Points:")
for p in transformed_points:
    print(f"[{p[0]}, {p[1]}]")

# Part 2: Visualize the effect on a unit square
def plot_transformation(T):
    # Create a figure with two subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Define the unit square
    square = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
    
    # Transform the square
    transformed_square = np.dot(square, T)
    
    # Plot the original square
    ax1.plot(square[:, 0], square[:, 1], 'b-')
    ax1.fill(square[:, 0], square[:, 1], 'lightblue', alpha=0.5)
    ax1.set_xlim(-0.5, 3.5)
    ax1.set_ylim(-0.5, 3.5)
    ax1.grid(True)
    ax1.set_title('Original Unit Square')
    ax1.set_aspect('equal')
    
    # Plot the transformed square
    ax2.plot(transformed_square[:, 0], transformed_square[:, 1], 'r-')
    ax2.fill(transformed_square[:, 0], transformed_square[:, 1], 'lightcoral', alpha=0.5)
    ax2.set_xlim(-0.5, 3.5)
    ax2.set_ylim(-0.5, 3.5)
    ax2.grid(True)
    ax2.set_title('Transformed Square')
    ax2.set_aspect('equal')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation
plot_transformation(T)

# Part 3: Find eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(T)

print("\nEigenvalues:")
for i, val in enumerate(eigenvalues):
    print(f"λ_{i+1} = {val:.4f}")

print("\nEigenvectors (as columns):")
print(eigenvectors)

# Part 4: Visualize eigenvectors and their transformations
def plot_eigenvectors(T, eigenvalues, eigenvectors):
    # Create a figure
    fig, ax = plt.subplots(figsize=(8, 8))
    
    # Set limits and grid
    ax.set_xlim(-3, 3)
    ax.set_ylim(-3, 3)
    ax.grid(True)
    ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    ax.axvline(x=0, color='k', linestyle='-', alpha=0.3)
    
    # Colors for the eigenvectors
    colors = ['blue', 'green']
    
    # Plot eigenvectors and their transformations
    for i in range(len(eigenvalues)):
        # Get the i-th eigenvector
        v = eigenvectors[:, i]
        
        # Scale eigenvector for better visualization
        v_scaled = v * 2
        
        # Plot the original eigenvector
        ax.arrow(0, 0, v_scaled[0], v_scaled[1], head_width=0.1, head_length=0.1, 
                 fc=colors[i], ec=colors[i], label=f"Eigenvector {i+1}")
        
        # Transform the eigenvector
        tv = np.dot(T, v)
        tv_scaled = tv * 2
        
        # Plot the transformed eigenvector
        ax.arrow(0, 0, tv_scaled[0], tv_scaled[1], head_width=0.1, head_length=0.1, 
                 fc='lightcoral', ec='lightcoral', linestyle='--')
        
        # Add text for eigenvalue
        ax.text(v_scaled[0]*1.1, v_scaled[1]*1.1, 
                f"λ_{i+1}={eigenvalues[i]:.2f}", 
                color=colors[i])
    
    # Add a unit circle and its transformation
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    
    # Stack coordinates to create points
    circle_points = np.column_stack((circle_x, circle_y))
    
    # Transform the circle
    transformed_circle = np.dot(circle_points, T)
    
    # Plot the unit circle
    ax.plot(circle_x, circle_y, 'k-', alpha=0.3)
    
    # Plot the transformed circle
    ax.plot(transformed_circle[:, 0], transformed_circle[:, 1], 'r-', alpha=0.3)
    
    ax.set_title('Eigenvectors and Their Transformations')
    ax.legend()
    plt.axis('equal')
    plt.tight_layout()
    plt.show()

# Visualize eigenvectors
plot_eigenvectors(T, eigenvalues, eigenvectors)

# Part 5: Visualize the transformation in 3D
def plot_transformation_3d():
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a grid of points
    x = np.linspace(-2, 2, 10)
    y = np.linspace(-2, 2, 10)
    X, Y = np.meshgrid(x, y)
    
    # Flatten the grid points
    points = np.column_stack((X.flatten(), Y.flatten()))
    
    # Transform the points
    transformed_points = np.dot(points, T)
    
    # Reshape the transformed coordinates
    X_transformed = transformed_points[:, 0].reshape(X.shape)
    Y_transformed = transformed_points[:, 1].reshape(Y.shape)
    
    # Create Z coordinates (all zeros for 2D transformation)
    Z = np.zeros(X.shape)
    Z_transformed = np.ones(X.shape)  # Offset for visualization
    
    # Plot the original grid
    ax.plot_surface(X, Y, Z, alpha=0.5, color='blue', label='Original')
    
    # Plot the transformed grid
    ax.plot_surface(X_transformed, Y_transformed, Z_transformed, alpha=0.5, color='red', label='Transformed')
    
    # Add lines connecting original and transformed points
    for i in range(len(points)):
        ax.plot([points[i, 0], transformed_points[i, 0]],
                [points[i, 1], transformed_points[i, 1]],
                [0, 1], 'k-', alpha=0.2)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z (Visualization Only)')
    ax.set_title('3D Visualization of 2D Linear Transformation')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation in 3D
plot_transformation_3d()

import numpy as np
import matplotlib.pyplot as plt

# Define the logistic map function
def logistic_map(r, x):
    return r * x * (1 - x)

# Parameters
n_iterations = 1000  # Total iterations
n_transient = 800    # Discard initial transient iterations
n_points = 200       # Number of points to plot per r value
r_values = np.linspace(2.5, 4.0, 1000)  # Range of r values
x0 = 0.1  # Initial condition (must be between 0 and 1)

# Arrays to store results
r_plot = []
x_plot = []

# Iterate over r values
for r in r_values:
    x = x0
    # Run iterations, discard transients
    for i in range(n_iterations):
        x = logistic_map(r, x)
        # Store the last n_points after transient
        if i >= n_transient:
            r_plot.append(r)
            x_plot.append(x)

# Create the bifurcation diagram
plt.figure(figsize=(10, 6))
plt.plot(r_plot, x_plot, ',k', alpha=0.25)  # Plot as small black dots
plt.xlabel('$r$')
plt.ylabel('$x$')
plt.title('Bifurcation Diagram of the Logistic Map')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Define the torus parametrization
def torus(u, v, R=2, r=1):
    """
    Parametrize a torus with major radius R and minor radius r
    u, v are the parameters in [0, 2π)
    """
    x = (R + r * np.cos(v)) * np.cos(u)
    y = (R + r * np.cos(v)) * np.sin(u)
    z = r * np.sin(v)
    return x, y, z

# Compute the metric coefficients for the torus
def torus_metric(u, v, R=2, r=1):
    """
    Compute the metric coefficients E, F, G for the torus at point (u, v)
    """
    E = (R + r * np.cos(v))**2  # g_11
    F = 0                       # g_12 = g_21
    G = r**2                    # g_22
    return E, F, G

# Compute the Christoffel symbols for the torus
def christoffel_symbols(u, v, R=2, r=1):
    """
    Compute the Christoffel symbols for the torus at point (u, v)
    Returns a 2x2x2 array where Gamma[i,j,k] = Gamma^i_{jk}
    """
    E, F, G = torus_metric(u, v, R, r)
    
    # Derivatives of metric coefficients
    dE_du = 0
    dE_dv = -2 * (R + r * np.cos(v)) * r * np.sin(v)
    dF_du = dF_dv = 0
    dG_du = 0
    dG_dv = 0
    
    # Initialize Christoffel symbols
    Gamma = np.zeros((2, 2, 2))
    
    # g^{11} = 1/E, g^{12} = g^{21} = 0, g^{22} = 1/G
    g_inv_11 = 1/E
    g_inv_22 = 1/G
    
    # Gamma^1_{11}
    Gamma[0, 0, 0] = 0.5 * g_inv_11 * dE_du
    
    # Gamma^1_{12} = Gamma^1_{21}
    Gamma[0, 0, 1] = Gamma[0, 1, 0] = 0.5 * g_inv_11 * dE_dv
    
    # Gamma^1_{22}
    Gamma[0, 1, 1] = -0.5 * g_inv_11 * dG_du
    
    # Gamma^2_{11}
    Gamma[1, 0, 0] = -0.5 * g_inv_22 * dE_dv
    
    # Gamma^2_{12} = Gamma^2_{21}
    Gamma[1, 0, 1] = Gamma[1, 1, 0] = 0.5 * g_inv_22 * dF_du
    
    # Gamma^2_{22}
    Gamma[1, 1, 1] = 0.5 * g_inv_22 * dG_dv
    
    return Gamma

# Geodesic equation as a first-order system
def geodesic_equation(t, y, R=2, r=1):
    """
    The geodesic equation as a first-order system
    y = [u, v, u_dot, v_dot]
    Returns [u_dot, v_dot, u_ddot, v_ddot]
    """
    u, v, u_dot, v_dot = y
    
    # Compute Christoffel symbols
    Gamma = christoffel_symbols(u, v, R, r)
    
    # Geodesic equations (second derivatives)
    u_ddot = -Gamma[0, 0, 0] * u_dot**2 - 2 * Gamma[0, 0, 1] * u_dot * v_dot - Gamma[0, 1, 1] * v_dot**2
    v_ddot = -Gamma[1, 0, 0] * u_dot**2 - 2 * Gamma[1, 0, 1] * u_dot * v_dot - Gamma[1, 1, 1] * v_dot**2
    
    return [u_dot, v_dot, u_ddot, v_ddot]

# Solve the geodesic equation with given initial conditions
def solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R=2, r=1):
    """
    Solve the geodesic equation with initial conditions
    u0, v0: initial position
    u_dot0, v_dot0: initial velocity
    t_span: time span [t_start, t_end]
    """
    y0 = [u0, v0, u_dot0, v_dot0]
    
    # Solve the ODE system
    sol = solve_ivp(
        lambda t, y: geodesic_equation(t, y, R, r),
        t_span,
        y0,
        method='RK45',
        rtol=1e-8,
        atol=1e-8,
        dense_output=True
    )
    
    return sol

# Plot the torus and geodesic
def plot_torus_and_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Plot the torus surface and the geodesic curve
    """
    # Create figure
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a mesh grid for the torus
    u = np.linspace(0, 2*np.pi, 100)
    v = np.linspace(0, 2*np.pi, 100)
    u_grid, v_grid = np.meshgrid(u, v)
    
    # Compute the torus coordinates
    x, y, z = torus(u_grid, v_grid, R, r)
    
    # Plot the torus surface with translucent appearance
    ax.plot_surface(x, y, z, cmap=cm.viridis, alpha=0.6, edgecolor='none')
    
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    
    # Convert the geodesic curve to Cartesian coordinates
    x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
    
    # Plot the geodesic curve
    ax.plot(x_geo, y_geo, z_geo, 'r-', linewidth=2, label='Geodesic')
    
    # Mark the start and end points
    ax.plot([x_geo[0]], [y_geo[0]], [z_geo[0]], 'go', markersize=8, label='Start')
    ax.plot([x_geo[-1]], [y_geo[-1]], [z_geo[-1]], 'bo', markersize=8, label='End')
    
    # Set labels and title
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Geodesic on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig, ax

# Analyze the length and curvature of the geodesic
def analyze_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Analyze properties of the geodesic such as length and curvature
    """
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    u_dot_sol, v_dot_sol = y_interp[2], y_interp[3]
    
    # Compute the metric coefficients along the geodesic
    E, F, G = [], [], []
    for i in range(len(t)):
        E_i, F_i, G_i = torus_metric(u_sol[i], v_sol[i], R, r)
        E.append(E_i)
        F.append(F_i)
        G.append(G_i)
    
    E = np.array(E)
    F = np.array(F)
    G = np.array(G)
    
    # Compute the speed along the geodesic
    speed = np.sqrt(E * u_dot_sol**2 + 2*F * u_dot_sol * v_dot_sol + G * v_dot_sol**2)
    
    # Compute the length of the geodesic
    dt = t[1] - t[0]
    length = np.sum(speed * dt)
    
    # Plot the speed along the geodesic
    plt.figure(figsize=(10, 6))
    plt.plot(t, speed)
    plt.title('Speed along the Geodesic')
    plt.xlabel('t')
    plt.ylabel('Speed')
    plt.grid(True)
    plt.show()
    
    print(f"Geodesic Length: {length:.6f}")
    
    return length

# Main execution code
if __name__ == "__main__":
    # Parameters for the torus
    R = 2.0  # Major radius
    r = 1.0  # Minor radius
    
    # Initial conditions for the geodesic
    u0 = 0.0
    v0 = 0.0
    u_dot0 = 1.0
    v_dot0 = 0.5
    
    # Time span for integration
    t_span = [0, 10]
    
    # Solve the geodesic equation
    sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
    
    # Plot the results
    plot_torus_and_geodesic(sol, R, r)
    
    # Analyze the geodesic
    geodesic_length = analyze_geodesic(sol, R, r)
    
    # Plot a few different geodesics
    plt.figure(figsize=(12, 10))
    
    initial_conditions = [
        (0.0, 0.0, 1.0, 0.0, 'r', 'Horizontal'),
        (0.0, 0.0, 0.0, 1.0, 'g', 'Vertical'),
        (0.0, 0.0, 1.0, 0.5, 'b', 'Diagonal'),
        (0.0, 0.0, 1.0, 2.0, 'purple', 'Steep Diagonal'),
    ]
    
    ax = plt.axes(projection='3d')
    
    # Plot the torus wireframe
    u = np.linspace(0, 2*np.pi, 30)
    v = np.linspace(0, 2*np.pi, 30)
    u_grid, v_grid = np.meshgrid(u, v)
    x, y, z = torus(u_grid, v_grid, R, r)
    ax.plot_surface(x, y, z, color='lightgray', alpha=0.2)
    
    # Compute and plot different geodesics
    for u0, v0, u_dot0, v_dot0, color, label in initial_conditions:
        sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
        
        # Extract the solution
        t = np.linspace(sol.t[0], sol.t[-1], 1000)
        y_interp = sol.sol(t)
        u_sol, v_sol = y_interp[0], y_interp[1]
        
        # Convert to Cartesian coordinates
        x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
        
        # Plot the geodesic
        ax.plot(x_geo, y_geo, z_geo, color=color, linewidth=2, label=label)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Different Geodesics on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Set random seed for reproducibility
np.random.seed(42)

# Parameters
n_samples = 10000  # Number of samples to generate
sample_sizes = [1, 2, 5, 10, 30, 50]  # Different sample sizes to demonstrate CLT
distributions = {
    'uniform': stats.uniform(0, 1),
    'exponential': stats.expon(scale=1),
    'poisson': stats.poisson(5)
}

# Function to generate samples and calculate means
def generate_sample_means(distribution, sample_size, n_samples):
    # Generate n_samples sets of sample_size samples and calculate their means
    sample_means = np.array([
        np.mean(distribution.rvs(size=sample_size)) 
        for _ in range(n_samples)
    ])
    return sample_means

# Setup the plots
fig, axes = plt.subplots(len(distributions), len(sample_sizes), figsize=(20, 12))
fig.suptitle('Central Limit Theorem Demonstration', fontsize=16)

# Create plots for each distribution and sample size
for i, (dist_name, distribution) in enumerate(distributions.items()):
    # Calculate the theoretical mean and std of the distribution
    if dist_name == 'uniform':
        true_mean = 0.5
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_mean = 1
        true_std = 1
    else:  # poisson
        true_mean = 5
        true_std = np.sqrt(5)
    
    for j, sample_size in enumerate(sample_sizes):
        # Generate sample means
        means = generate_sample_means(distribution, sample_size, n_samples)
        
        # Calculate expected standard deviation of the sampling distribution
        expected_std = true_std / np.sqrt(sample_size)
        
        # Plot histogram
        ax = axes[i, j]
        ax.hist(means, bins=50, density=True, alpha=0.6, color='skyblue')
        
        # Plot the normal distribution curve for comparison
        x = np.linspace(min(means), max(means), 1000)
        normal_pdf = stats.norm.pdf(x, true_mean, expected_std)
        ax.plot(x, normal_pdf, 'r-', linewidth=2)
        
        # Add theoretical mean line
        ax.axvline(true_mean, color='green', linestyle='--', linewidth=1.5)
        
        # Set title and labels
        if i == 0:
            ax.set_title(f'Sample Size = {sample_size}')
        if j == 0:
            ax.set_ylabel(f'{dist_name.capitalize()}')
        
        # Add some stats
        sample_mean = np.mean(means)
        sample_std = np.std(means)
        ax.text(0.05, 0.95, f'Mean: {sample_mean:.4f}\nSTD: {sample_std:.4f}', 
                transform=ax.transAxes, fontsize=8, 
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))

plt.tight_layout()
plt.subplots_adjust(top=0.9)

# Create additional visualization for convergence of variance
plt.figure(figsize=(10, 6))
colors = ['blue', 'red', 'green']
markers = ['o', 's', '^']

for i, (dist_name, distribution) in enumerate(distributions.items()):
    if dist_name == 'uniform':
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_std = 1
    else:  # poisson
        true_std = np.sqrt(5)
    
    std_ratios = []
    for sample_size in range(1, 51):
        means = generate_sample_means(distribution, sample_size, n_samples)
        sample_std = np.std(means)
        expected_std = true_std / np.sqrt(sample_size)
        std_ratio = sample_std / expected_std
        std_ratios.append(std_ratio)
    
    plt.plot(range(1, 51), std_ratios, color=colors[i], marker=markers[i], 
             linestyle='-', label=f'{dist_name.capitalize()}', markersize=4)

plt.axhline(1, color='black', linestyle='--')
plt.xlabel('Sample Size')
plt.ylabel('Ratio of Sample STD to Expected STD')
plt.title('Convergence of Sampling Distribution Standard Deviation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Let's also visualize the original distributions
plt.figure(figsize=(15, 5))
x_values = {
    'uniform': np.linspace(0, 1, 1000),
    'exponential': np.linspace(0, 5, 1000),
    'poisson': np.arange(0, 20)
}

for i, (dist_name, distribution) in enumerate(distributions.items()):
    plt.subplot(1, 3, i+1)
    
    if dist_name == 'poisson':
        # For discrete distribution
        x = x_values[dist_name]
        plt.bar(x, distribution.pmf(x), alpha=0.6, color=colors[i])
    else:
        # For continuous distributions
        x = x_values[dist_name]
        plt.plot(x, distribution.pdf(x), color=colors[i], linewidth=2)
        plt.fill_between(x, distribution.pdf(x), alpha=0.3, color=colors[i])
    
    plt.title(f'{dist_name.capitalize()} Distribution')
    plt.xlabel('Value')
    plt.ylabel('Density' if dist_name != 'poisson' else 'Probability')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()