具象直观地解释这一段落内容。

import pandas as pd
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm

# 创建数据
data = {
    '面积': [80, 95, 110, 120, 135, 145, 160],
    '卧室数': [2, 3, 3, 4, 4, 4, 5],
    '房龄': [5, 3, 8, 10, 2, 6, 12],
    '房价': [320, 420, 480, 520, 580, 620, 680]
}
df = pd.DataFrame(data)

# 建立模型
X = df[['面积', '卧室数', '房龄']]
y = df['房价']
X = sm.add_constant(X)  # 添加截距项

model = sm.OLS(y, X)
results = model.fit()

print(results.summary())

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures

# 生长数据：时间(周) vs 高度(cm)
时间 = np.array([1, 2, 3, 4, 5, 6, 7, 8])
高度 = np.array([2, 5, 8, 10, 11, 11.5, 11.8, 12])

# 使用2阶多项式
poly = PolynomialFeatures(degree=2)
X_poly = poly.fit_transform(时间.reshape(-1, 1))

model = LinearRegression()
model.fit(X_poly, 高度)

# 预测和绘图
X_fit = np.linspace(0, 9, 100).reshape(-1, 1)
X_fit_poly = poly.transform(X_fit)
y_fit = model.predict(X_fit_poly)

plt.scatter(时间, 高度, color='blue', label='实际数据')
plt.plot(X_fit, y_fit, color='red', label='多项式拟合')
plt.xlabel('时间（周）')
plt.ylabel('高度（cm）')
plt.legend()
plt.show()

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix

# 信用数据
data = {
    '收入': [5000, 8000, 3000, 6000, 4000, 7000, 5500, 3500],
    '负债比': [0.2, 0.1, 0.6, 0.3, 0.5, 0.2, 0.4, 0.7],
    '违约': [0, 0, 1, 0, 1, 0, 0, 1]  # 0:不违约, 1:违约
}
df = pd.DataFrame(data)

# 建立逻辑回归模型
X = df[['收入', '负债比']]
y = df['违约']

model = LogisticRegression()
model.fit(X, y)

# 预测
y_pred = model.predict(X)
y_prob = model.predict_proba(X)[:, 1]  # 预测为1类的概率

print(f"准确率: {accuracy_score(y, y_pred):.3f}")
print("混淆矩阵:")
print(confusion_matrix(y, y_pred))

# 查看系数
print(f"系数: {model.coef_}")
print(f"截距: {model.intercept_}")

from statsmodels.tsa.stattools import adfuller

result = adfuller(data)
p_value = result[1]
is_stationary = p_value < 0.05

import matplotlib.pyplot as plt
plt.plot(data)
plt.show()

from statsmodels.tsa.arima.model import ARIMA

model = ARIMA(data, order=(1,1,1))
results = model.fit()

# 预测未来7天
forecast = results.get_forecast(steps=7)
mean = forecast.predicted_mean          # 点预测
conf_int = forecast.conf_int()          # 置信区间

import numpy as np

class AHP:
    def __init__(self, comparison_matrix):
        self.A = np.array(comparison_matrix)
        self.n = self.A.shape[0]
        self.RI = {1:0, 2:0, 3:0.52, 4:0.89, 5:1.12, 6:1.26, 7:1.36, 8:1.41, 9:1.46}

    def calculate_weights(self):
        # 几何平均法计算权重
        geometric_means = np.prod(self.A, axis=1)  (1/self.n)
        weights = geometric_means / np.sum(geometric_means)
        return weights

    def consistency_check(self, weights):
        # 计算最大特征值
        AW = np.dot(self.A, weights)
        lambda_max = np.mean(AW / weights)

        # 一致性指标
        CI = (lambda_max - self.n) / (self.n - 1)
        CR = CI / self.RI[self.n]

        return lambda_max, CI, CR

    def analyze(self):
        weights = self.calculate_weights()
        lambda_max, CI, CR = self.consistency_check(weights)

        print("判断矩阵:")
        print(self.A)
        print(f"\n权重向量: {weights}")
        print(f"最大特征值: {lambda_max:.4f}")
        print(f"一致性指标 CI: {CI:.4f}")
        print(f"一致性比率 CR: {CR:.4f}")

        if CR < 0.1:
            print("✓ 一致性检验通过")
        else:
            print("✗ 一致性检验未通过，需要调整判断矩阵")

        return weights, CR

# 判断矩阵
comparison_matrix = [
    [1, 1/3, 3, 2],      # 价格 vs 其他
    [3, 1, 5, 4],        # 性能 vs 其他  
    [1/3, 1/5, 1, 1/2],  # 外观 vs 其他
    [1/2, 1/4, 2, 1]     # 品牌 vs 其他
]

# 执行分析
ahp = AHP(comparison_matrix)
weights, cr = ahp.analyze()

判断矩阵:
[[1.         0.33333333 3.         2.        ]
 [3.         1.         5.         4.        ]
 [0.33333333 0.2        1.         0.5       ]
 [0.5        0.25       2.         1.        ]]

权重向量: [0.245 0.476 0.086 0.193]
最大特征值: 4.0512
一致性指标 CI: 0.0171
一致性比率 CR: 0.0190
✓ 一致性检验通过

import numpy as np
import pandas as pd

def entropy_weight(data):
    """
    计算指标的熵权
    :param data: DataFrame，行为样本，列为指标。所有指标必须为正向指标且非负。
    :return: 权重向量
    """
    # 步骤1：数据标准化 (比重法)
    data = data / data.sum(axis=0)

    # 步骤2：计算信息熵
    data = data.replace(0, 1e-10) # 避免log(0)
    m = data.shape[0]
    k = 1.0 / np.log(m)
    e = -k * (data * np.log(data)).sum(axis=0)

    # 步骤3 & 4：计算权重
    d = 1 - e
    w = d / d.sum()
    return w.values

# 案例应用
data_df = pd.DataFrame({
    '成绩': [90, 80, 70, 60],
    '纪律得分': [1, 2, 3, 0],
    '班干部数量': [1, 2, 3, 0]
})

print("原始数据矩阵：")
print(data_df)

weights = entropy_weight(data_df)
print("\n各指标熵权：")
for col, w in zip(data_df.columns, weights):
    print(f"- {col}: {w:.4f} ({w*100:.2f}%)")

import pulp
import numpy as np

def linear_programming_example():
    """
    生产计划优化问题实例
    目标：在资源约束下最大化利润
    问题描述：
      工厂生产两种产品A和B
      - 产品A：利润80元/件，耗时2小时，耗料4单位
      - 产品B：利润60元/件，耗时3小时，耗料2单位
      资源限制：
      - 每日可用工时：100小时
      - 每日可用原料：80单位
      - 市场需求：产品A不超过30件，产品B不超过20件
    """
    # 创建优化问题，指定问题名称和优化方向（最大化）
    prob = pulp.LpProblem("Production_Planning", pulp.LpMaximize)
    
    # 定义决策变量
    # x1: 产品A的日产量，下限0，上限30，整数约束
    # x2: 产品B的日产量，下限0，上限20，整数约束
    x1 = pulp.LpVariable("Product_A", lowBound=0, upBound=30, cat='Integer')
    x2 = pulp.LpVariable("Product_B", lowBound=0, upBound=20, cat='Integer')
    
    # 设置目标函数：最大化总利润
    prob += 80*x1 + 60*x2, "Total_Profit"
    
    # 添加约束条件
    prob += 2*x1 + 3*x2 <= 100, "Labor_Constraint"    # 工时约束
    prob += 4*x1 + 2*x2 <= 80, "Material_Constraint"  # 材料约束
    
    # 求解问题
    prob.solve()
    
    # 输出结果
    print("线性规划求解结果：")
    print(f"求解状态: {pulp.LpStatus[prob.status]}")
    print(f"最优解: 生产产品A {x1.varValue} 件，产品B {x2.varValue} 件")
    print(f"最大利润: {pulp.value(prob.objective)} 元")
    print(f"工时资源使用: {2*x1.varValue + 3*x2.varValue} / 100 小时")
    print(f"材料资源使用: {4*x1.varValue + 2*x2.varValue} / 80 单位")
    
    return x1.varValue, x2.varValue

# 执行示例
if __name__ == "__main__":
    optimal_A, optimal_B = linear_programming_example()

import numpy as np
import random
from typing import List, Tuple, Callable

class GeneticAlgorithm:
    """
    遗传算法实现类
    用于解决函数优化问题
    """
    
    def __init__(self, population_size: int, mutation_rate: float, 
                 crossover_rate: float, generations: int):
        self.population_size = population_size
        self.mutation_rate = mutation_rate
        self.crossover_rate = crossover_rate
        self.generations = generations
        
    def initialize_population(self, bounds: List[Tuple[float, float]]) -> List[List[float]]:
        """初始化种群，在给定边界内随机生成个体"""
        population = []
        for _ in range(self.population_size):
            individual = [random.uniform(low, high) for low, high in bounds]
            population.append(individual)
        return population
    
    def evaluate_fitness(self, population: List[List[float]], 
                        fitness_func: Callable) -> List[float]:
        """评估种群中每个个体的适应度"""
        return [fitness_func(ind) for ind in population]
    
    def selection(self, population: List[List[float]], 
                 fitnesses: List[float]) -> List[List[float]]:
        """锦标赛选择操作"""
        selected = []
        for _ in range(len(population)):
            # 随机选择3个个体进行锦标赛
            tournament_indices = random.sample(range(len(population)), 3)
            tournament_fitness = [fitnesses[i] for i in tournament_indices]
            winner_idx = tournament_indices[np.argmax(tournament_fitness)]
            selected.append(population[winner_idx])
        return selected
    
    def crossover(self, parent1: List[float], parent2: List[float]) -> Tuple[List[float]]:
        """算术交叉操作"""
        if random.random() < self.crossover_rate:
            alpha = random.random()
            child1 = [alpha * p1 + (1 - alpha) * p2 for p1, p2 in zip(parent1, parent2)]
            child2 = [alpha * p2 + (1 - alpha) * p1 for p1, p2 in zip(parent1, parent2)]
            return child1, child2
        return parent1, parent2
    
    def mutation(self, individual: List[float], bounds: List[Tuple[float, float]]) -> List[float]:
        """高斯变异操作"""
        mutated = individual.copy()
        for i in range(len(mutated)):
            if random.random() < self.mutation_rate:
                low, high = bounds[i]
                # 在当前值基础上添加高斯噪声
                mutated[i] += random.gauss(0, (high - low) * 0.1)
                # 确保变异后的值仍在边界内
                mutated[i] = max(low, min(high, mutated[i]))
        return mutated
    
    def run(self, fitness_func: Callable, bounds: List[Tuple[float, float]]) -> dict:
        """运行遗传算法"""
        # 初始化
        population = self.initialize_population(bounds)
        best_individual = None
        best_fitness = -float('inf')
        fitness_history = []
        
        for generation in range(self.generations):
            # 评估适应度
            fitnesses = self.evaluate_fitness(population, fitness_func)
            
            # 更新历史最佳
            current_best_fitness = max(fitnesses)
            if current_best_fitness > best_fitness:
                best_fitness = current_best_fitness
                best_individual = population[np.argmax(fitnesses)]
            
            fitness_history.append(best_fitness)
            
            # 遗传操作
            selected = self.selection(population, fitnesses)
            
            # 交叉
            new_population = []
            for i in range(0, len(selected), 2):
                if i + 1 < len(selected):
                    child1, child2 = self.crossover(selected[i], selected[i+1])
                    new_population.extend([child1, child2])
                else:
                    new_population.append(selected[i])
            
            # 变异
            population = [self.mutation(ind, bounds) for ind in new_population]
            
            # 输出进度
            if generation % 50 == 0:
                print(f"Generation {generation}: Best Fitness = {best_fitness:.6f}")
        
        return {
            'best_individual': best_individual,
            'best_fitness': best_fitness,
            'fitness_history': fitness_history
        }

# 使用示例：优化Rastrigin函数
def rastrigin_function(x: List[float]) -> float:
    """Rastrigin测试函数，具有多个局部极小值"""
    A = 10
    n = len(x)
    return - (A * n + sum(xi2 - A * np.cos(2 * np.pi * xi) for xi in x))

# 执行遗传算法优化
if __name__ == "__main__":
    # 设置算法参数
    ga = GeneticAlgorithm(
        population_size=50,
        mutation_rate=0.1,
        crossover_rate=0.8,
        generations=200
    )
    
    # 定义搜索空间（2维问题）
    bounds = [(-5.12, 5.12), (-5.12, 5.12)]
    
    # 运行优化
    result = ga.run(rastrigin_function, bounds)
    
    print("\n遗传算法优化结果：")
    print(f"最优解: {result['best_individual']}")
    print(f"最优值: {-result['best_fitness']:.6f}")  # 注意我们求的是负函数值的最小化

def knapsack_dp(weights: List[int], values: List[int], capacity: int) -> Tuple[int, List[int]]:
    """
    0-1背包问题动态规划求解
    参数:
        weights: 物品重量列表
        values: 物品价值列表
        capacity: 背包容量
    返回:
        最大总价值和选择的物品索引列表
    """
    n = len(weights)
    
    # 创建DP表，dp[i][w]表示前i个物品在容量w下的最大价值
    dp = [[0] * (capacity + 1) for _ in range(n + 1)]
    
    # 填充DP表
    for i in range(1, n + 1):
        for w in range(1, capacity + 1):
            if weights[i-1] <= w:
                # 可以选择放入或不放入当前物品
                dp[i][w] = max(dp[i-1][w], 
                              values[i-1] + dp[i-1][w - weights[i-1]])
            else:
                # 无法放入当前物品
                dp[i][w] = dp[i-1][w]
    
    # 回溯找出选择的物品
    selected_items = []
    w = capacity
    for i in range(n, 0, -1):
        if dp[i][w] != dp[i-1][w]:
            selected_items.append(i-1)
            w -= weights[i-1]
    
    # 反转列表使物品按原始顺序排列
    selected_items.reverse()
    
    return dp[n][capacity], selected_items

def investment_planning_dp():
    """
    投资规划问题动态规划求解
    问题描述：
      有100万元资金，可以投资到3个项目
      每个项目有不同的投资额和收益
      目标：最大化总收益
    """
    # 投资选项：每个元组表示(投资额, 收益)
    investments = [
        (0, 0),      # 不投资
        (20, 15),    # 项目1：投资20万，收益15万
        (30, 25),    # 项目2：投资30万，收益25万  
        (50, 40),    # 项目3：投资50万，收益40万
        (70, 55)     # 项目4：投资70万，收益55万
    ]
    
    total_capital = 100  # 总资金100万元
    n = len(investments)
    
    # DP表：dp[i][j]表示考虑前i个项目，使用j万元资金的最大收益
    dp = [[0] * (total_capital + 1) for _ in range(n)]
    
    # 记录决策路径
    decision = [[0] * (total_capital + 1) for _ in range(n)]
    
    # 动态规划求解
    for i in range(1, n):
        investment, profit = investments[i]
        for j in range(total_capital + 1):
            # 不投资当前项目
            dp[i][j] = dp[i-1][j]
            decision[i][j] = 0
            
            # 如果资金足够，考虑投资当前项目
            if j >= investment:
                if dp[i-1][j - investment] + profit > dp[i][j]:
                    dp[i][j] = dp[i-1][j - investment] + profit
                    decision[i][j] = i
    
    # 回溯找出最优投资方案
    current_capital = total_capital
    investment_plan = []
    for i in range(n-1, 0, -1):
        if decision[i][current_capital] != 0:
            investment_plan.append(decision[i][current_capital])
            current_capital -= investments[decision[i][current_capital]][0]
    
    investment_plan.reverse()
    
    return dp[n-1][total_capital], investment_plan

# 执行示例
if __name__ == "__main__":
    print("动态规划算法示例")
    print("=" * 50)
    
    # 0-1背包问题示例
    weights = [2, 3, 4, 5, 9]
    values = [3, 4, 5, 8, 10]
    capacity = 20
    
    max_value, selected = knapsack_dp(weights, values, capacity)
    print("0-1背包问题求解结果：")
    print(f"物品重量: {weights}")
    print(f"物品价值: {values}")
    print(f"背包容量: {capacity}")
    print(f"最大价值: {max_value}")
    print(f"选择的物品索引: {selected}")
    print(f"选择物品的重量: {[weights[i] for i in selected]}")
    print(f"总重量: {sum(weights[i] for i in selected)}")
    
    print("\n" + "=" * 50)
    
    # 投资规划问题示例
    max_profit, plan = investment_planning_dp()
    print("投资规划问题求解结果：")
    print(f"可用资金: 100万元")
    print(f"最大收益: {max_profit}万元")
    print(f"投资方案: {plan}")
    
    # 计算总投资额
    investments = [(0,0), (20,15), (30,25), (50,40), (70,55)]
    total_investment = sum(investments[i][0] for i in plan)
    print(f"总投资额: {total_investment}万元")
    print(f"资金利用率: {total_investment}/100 = {total_investment}%")

import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

def kmeans_analysis(X, max_k=10):
    """
    K-Means聚类完整分析
    """
    # 数据标准化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # 肘部法则确定最优k值
    sse = []
    k_range = range(1, max_k + 1)
    
    for k in k_range:
        kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
        kmeans.fit(X_scaled)
        sse.append(kmeans.inertia_)
    
    # 绘制肘部图
    plt.figure(figsize=(10, 6))
    plt.plot(k_range, sse, 'bo-')
    plt.xlabel('Number of Clusters (k)')
    plt.ylabel('Sum of Squared Errors (SSE)')
    plt.title('Elbow Method for Optimal k')
    plt.grid(True)
    plt.show()
    
    # 计算SSE下降率确定肘点
    sse_reduction = []
    for i in range(1, len(sse) - 1):
        reduction = (sse[i-1] - sse[i]) / (sse[i] - sse[i+1])
        sse_reduction.append(reduction)
    
    optimal_k = np.argmax(sse_reduction) + 2  # +2因为从k=2开始计算
    
    print(f"根据肘部法则，推荐簇数: {optimal_k}")
    
    # 使用最优k值进行聚类
    kmeans_optimal = KMeans(n_clusters=optimal_k, random_state=42)
    labels = kmeans_optimal.fit_predict(X_scaled)
    centers = kmeans_optimal.cluster_centers_
    
    return labels, centers, kmeans_optimal

# 使用示例
# 假设 X 是 n×d 的数据矩阵
# labels, centers, model = kmeans_analysis(X)

from sklearn.cluster import AgglomerativeClustering
import scipy.cluster.hierarchy as sch
import matplotlib.pyplot as plt

def hierarchical_clustering_analysis(X, method='ward'):
    """
    层次聚类完整分析
    """
    # 数据标准化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # 绘制树状图
    plt.figure(figsize=(12, 8))
    dendrogram = sch.dendrogram(
        sch.linkage(X_scaled, method=method),
        truncate_mode='lastp',
        p=30,
        show_leaf_counts=True,
        leaf_rotation=90,
        leaf_font_size=10
    )
    plt.title(f'Hierarchical Clustering Dendrogram ({method} linkage)')
    plt.xlabel('Sample Index')
    plt.ylabel('Distance')
    plt.show()
    
    # 从树状图确定最佳簇数（示例）
    # 在实际应用中，应根据树状图的垂直距离跳跃确定切割高度
    optimal_k = 3  # 这需要根据具体树状图确定
    
    # 执行层次聚类
    agg_clustering = AgglomerativeClustering(
        n_clusters=optimal_k, 
        linkage=method
    )
    labels = agg_clustering.fit_predict(X_scaled)
    
    return labels, dendrogram

# 使用示例
# labels, dendrogram = hierarchical_clustering_analysis(X, method='ward')

from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors
import numpy as np

def dbscan_analysis(X, min_samples=None):
    """
    DBSCAN聚类完整分析
    """
    # 数据标准化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    
    # 自动确定min_samples
    if min_samples is None:
        min_samples = 2 * X_scaled.shape[1]
    
    # 寻找最优eps参数
    neighbors = NearestNeighbors(n_neighbors=min_samples)
    neighbors_fit = neighbors.fit(X_scaled)
    distances, indices = neighbors_fit.kneighbors(X_scaled)
    
    k_distances = np.sort(distances[:, min_samples-1])
    
    # 绘制k-距离图
    plt.figure(figsize=(10, 6))
    plt.plot(k_distances)
    plt.xlabel('Data Points Sorted by Distance')
    plt.ylabel(f'{min_samples}-NN Distance')
    plt.title('K-Distance Graph for Optimal Eps')
    plt.grid(True)
    
    # 自动选择eps（使用第95百分位数）
    optimal_eps = np.percentile(k_distances, 95)
    plt.axhline(y=optimal_eps, color='r', linestyle='--', 
                label=f'Recommended eps: {optimal_eps:.2f}')
    plt.legend()
    plt.show()
    
    # 执行DBSCAN聚类
    dbscan = DBSCAN(eps=optimal_eps, min_samples=min_samples)
    labels = dbscan.fit_predict(X_scaled)
    
    # 统计结果
    n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
    n_noise = list(labels).count(-1)
    
    print(f"DBSCAN聚类结果:")
    print(f"簇数量: {n_clusters}")
    print(f"噪声点数量: {n_noise}")
    print(f"噪声点比例: {n_noise/len(X_scaled):.2%}")
    
    return labels, dbscan

# 使用示例
# labels, dbscan_model = dbscan_analysis(X)

from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score

def evaluate_clustering(X, labels):
    """
    综合评估聚类结果
    """
    if len(set(labels)) < 2:
        print("警告：簇数量少于2，无法计算评估指标")
        return None
    
    metrics = {}
    
    # 轮廓系数
    metrics['silhouette'] = silhouette_score(X, labels)
    
    # Calinski-Harabasz指数
    metrics['calinski_harabasz'] = calinski_harabasz_score(X, labels)
    
    # Davies-Bouldin指数
    metrics['davies_bouldin'] = davies_bouldin_score(X, labels)
    
    # 打印结果
    print("=== 聚类质量评估 ===")
    print(f"轮廓系数: {metrics['silhouette']:.4f}")
    print(f"Calinski-Harabasz指数: {metrics['calinski_harabasz']:.2f}")
    print(f"Davies-Bouldin指数: {metrics['davies_bouldin']:.4f}")
    
    # 轮廓系数解释
    silhouette = metrics['silhouette']
    if silhouette > 0.7:
        interpretation = "强聚类结构"
    elif silhouette > 0.5:
        interpretation = "合理的聚类结构"
    elif silhouette > 0.25:
        interpretation = "弱的聚类结构"
    else:
        interpretation = "没有实质的聚类结构"
    
    print(f"轮廓系数解释: {interpretation}")
    
    return metrics

# 使用示例
# metrics = evaluate_clustering(X_scaled, labels)

import numpy as np
import matplotlib.pyplot as plt

def monte_carlo_pi(num_samples):
    """使用蒙特卡洛方法计算圆周率π"""
    np.random.seed(42)
    
    # 在正方形区域内随机生成点
    x = np.random.uniform(-1, 1, num_samples)
    y = np.random.uniform(-1, 1, num_samples)
    
    # 计算每个点到原点的距离
    distance_squared = x2 + y2
    
    # 判断点是否在圆内
    inside_circle = distance_squared <= 1
    
    # 统计圆内的点数
    points_inside = np.sum(inside_circle)
    
    # 计算π的估计值
    pi_estimate = 4 * points_inside / num_samples
    
    return pi_estimate, x, y, inside_circle

# 使用不同样本数量进行测试
sample_sizes = [1000, 10000, 100000, 1000000]
results = {}

for n in sample_sizes:
    pi_est, x, y, inside = monte_carlo_pi(n)
    results[n] = pi_est
    error = abs(pi_est - np.pi)
    print(f"样本数: {n:8d} | π估计值: {pi_est:.6f} | 误差: {error:.6f}")

# 可视化最后一种情况
plt.figure(figsize=(8, 8))
plt.scatter(x[inside], y[inside], color='blue', s=1, alpha=0.6, label='圆内点')
plt.scatter(x[~inside], y[~inside], color='red', s=1, alpha=0.6, label='圆外点')
plt.title(f'蒙特卡洛模拟计算 π (样本数: {n}, 估计值: {pi_est:.4f})')
plt.axis('equal')
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def geometric_brownian_motion(S0, mu, sigma, T, dt, num_paths):
    """模拟几何布朗运动路径"""
    np.random.seed(42)
    
    num_steps = int(T / dt)
    
    # 初始化价格数组
    S = np.zeros((num_steps + 1, num_paths))
    S[0] = S0
    
    # 模拟路径
    for t in range(1, num_steps + 1):
        # 生成标准正态分布的随机数
        Z = np.random.standard_normal(num_paths)
        # 根据几何布朗运动公式更新价格
        S[t] = S[t-1] * np.exp((mu - 0.5 * sigma2) * dt + sigma * np.sqrt(dt) * Z)
    
    return S

# 参数设置
S0 = 100.0       # 初始股价
mu = 0.05        # 预期年化回报率 (5%)
sigma = 0.2      # 年化波动率 (20%)
T = 2.0          # 时间长度 (2年)
dt = 1/252       # 时间步长 (交易日)
num_paths = 10   # 模拟路径数量

# 运行模拟
S = geometric_brownian_motion(S0, mu, sigma, T, dt, num_paths)

# 时间轴
time = np.linspace(0, T, len(S))

# 绘制结果
plt.figure(figsize=(12, 6))

# 绘制所有路径
for i in range(num_paths):
    plt.plot(time, S[:, i], linewidth=1, alpha=0.7)

# 绘制均值路径
mean_path = np.mean(S, axis=1)
plt.plot(time, mean_path, 'k--', linewidth=2, label='平均路径')

plt.title('几何布朗运动模拟 - 股票价格路径')
plt.xlabel('时间 (年)')
plt.ylabel('股票价格')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

# 统计分析
final_prices = S[-1, :]
print("模拟结果统计:")
print(f"最终价格均值: {np.mean(final_prices):.2f}")
print(f"最终价格标准差: {np.std(final_prices):.2f}")
print(f"理论期望值: {S0 * np.exp(mu * T):.2f}")
print(f"最低最终价格: {np.min(final_prices):.2f}")
print(f"最高最终价格: {np.max(final_prices):.2f}")

def european_call_option_price(S0, K, T, r, sigma, num_simulations):
    """使用蒙特卡洛方法计算欧式看涨期权价格"""
    np.random.seed(42)
    
    # 模拟最终股价
    Z = np.random.standard_normal(num_simulations)
    ST = S0 * np.exp((r - 0.5 * sigma2) * T + sigma * np.sqrt(T) * Z)
    
    # 计算期权收益
    payoffs = np.maximum(ST - K, 0)
    
    # 计算现值
    option_price = np.exp(-r * T) * np.mean(payoffs)
    
    # 计算标准误
    std_error = np.std(payoffs) / np.sqrt(num_simulations)
    
    return option_price, std_error, payoffs

# 期权参数
S0 = 100    # 当前股价
K = 105     # 行权价
T = 1.0     # 到期时间（年）
r = 0.03    # 无风险利率
sigma = 0.2 # 波动率
num_simulations = 100000

# 计算期权价格
price, std_error, payoffs = european_call_option_price(S0, K, T, r, sigma, num_simulations)

print(f"欧式看涨期权定价结果:")
print(f"期权价格: {price:.4f}")
print(f"标准误差: {std_error:.6f}")
print(f"95% 置信区间: [{price - 1.96*std_error:.4f}, {price + 1.96*std_error:.4f}]")

# 可视化收益分布
plt.figure(figsize=(10, 6))
plt.hist(payoffs, bins=50, alpha=0.7, edgecolor='black')
plt.axvline(price * np.exp(r * T), color='red', linestyle='--', label='平均收益')
plt.title('期权收益分布')
plt.xlabel('收益')
plt.ylabel('频数')
plt.legend()
plt.show()