Python NumPy生成随机数据并进行统计分析详解

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子以确保结果可重现
np.random.seed(42)

# 生成均匀分布的随机数
uniform_data = np.random.uniform(0, 1, 1000)
print(f"均匀分布随机数样本: {uniform_data[:5]}")

# 生成正态分布的随机数
normal_data = np.random.normal(0, 1, 1000)
print(f"正态分布随机数样本: {normal_data[:5]}")

# 生成整数随机数
integer_data = np.random.randint(1, 100, 1000)
print(f"整数随机数样本: {integer_data[:5]}")

def simulate_dice_rolls(num_rolls):
    """
    模拟掷骰子实验
    
    参数:
    num_rolls: 掷骰子的次数
    
    返回:
    包含每次掷骰子结果的数组
    """
    # 生成1到6之间的随机整数
    dice_rolls = np.random.randint(1, 7, num_rolls)
    return dice_rolls

# 模拟1000次掷骰子
dice_results = simulate_dice_rolls(1000)

# 统计每个点数出现的频率
unique, counts = np.unique(dice_results, return_counts=true)
probability = counts / len(dice_results)

print("掷骰子结果统计:")
for i in range(len(unique)):
    print(f"点数 {unique[i]}: 出现 {counts[i]} 次, 概率 {probability[i]:.3f}")

# 计算基本统计量
mean_value = np.mean(dice_results)
std_dev = np.std(dice_results)
median_value = np.median(dice_results)

print(f"\n统计摘要:")
print(f"平均值: {mean_value:.2f}")
print(f"标准差: {std_dev:.2f}")
print(f"中位数: {median_value:.2f}")

def simulate_stock_prices(initial_price, mu, sigma, days, simulations):
    """
    使用几何布朗运动模拟股票价格
    
    参数:
    initial_price: 初始股价
    mu: 预期收益率
    sigma: 波动率
    days: 模拟天数
    simulations: 模拟次数
    
    返回:
    股价路径矩阵
    """
    # 生成时间步长
    dt = 1/252  # 假设一年有252个交易日
    
    # 生成随机数
    random_shocks = np.random.normal(0, 1, (days, simulations))
    
    # 初始化股价矩阵
    price_paths = np.zeros((days, simulations))
    price_paths[0] = initial_price
    
    # 计算每日收益率
    for t in range(1, days):
        drift = (mu - 0.5 * sigma**2) * dt
        diffusion = sigma * np.sqrt(dt) * random_shocks[t]
        price_paths[t] = price_paths[t-1] * np.exp(drift + diffusion)
    
    return price_paths

# 模拟参数
initial_price = 100  # 初始股价
mu = 0.1  # 年化预期收益率10%
sigma = 0.2  # 年化波动率20%
days = 252  # 一年的交易日
simulations = 1000  # 进行1000次模拟

# 执行模拟
stock_prices = simulate_stock_prices(initial_price, mu, sigma, days, simulations)

# 计算最终价格的统计信息
final_prices = stock_prices[-1]
mean_final_price = np.mean(final_prices)
std_final_price = np.std(final_prices)
min_final_price = np.min(final_prices)
max_final_price = np.max(final_prices)

print("股票价格模拟结果:")
print(f"初始价格: ${initial_price:.2f}")
print(f"平均最终价格: ${mean_final_price:.2f}")
print(f"价格标准差: ${std_final_price:.2f}")
print(f"最低最终价格: ${min_final_price:.2f}")
print(f"最高最终价格: ${max_final_price:.2f}")

# 计算价格变化百分比
price_changes = (final_prices - initial_price) / initial_price * 100
positive_changes = np.sum(price_changes > 0)
negative_changes = np.sum(price_changes < 0)

print(f"\n价格变化统计:")
print(f"价格上涨的模拟次数: {positive_changes} ({positive_changes/simulations*100:.1f}%)")
print(f"价格下跌的模拟次数: {negative_changes} ({negative_changes/simulations*100:.1f}%)")

def generate_temperature_data(days=365, base_temp=20, seasonal_variation=15, noise_level=3):
    """
    生成模拟的温度数据
    
    参数:
    days: 天数
    base_temp: 基础温度
    seasonal_variation: 季节性变化幅度
    noise_level: 随机噪声水平
    
    返回:
    温度数据数组
    """
    # 创建时间序列
    time = np.arange(days)
    
    # 添加季节性变化（使用正弦函数）
    seasonal_component = seasonal_variation * np.sin(2 * np.pi * time / 365)
    
    # 添加随机噪声
    noise = np.random.normal(0, noise_level, days)
    
    # 生成温度数据
    temperatures = base_temp + seasonal_component + noise
    
    return temperatures

# 生成一年的温度数据
temperature_data = generate_temperature_data()

# 基本统计分析
temp_mean = np.mean(temperature_data)
temp_std = np.std(temperature_data)
temp_min = np.min(temperature_data)
temp_max = np.max(temperature_data)
temp_median = np.median(temperature_data)

print("温度数据分析结果:")
print(f"平均温度: {temp_mean:.2f}°c")
print(f"温度标准差: {temp_std:.2f}°c")
print(f"最低温度: {temp_min:.2f}°c")
print(f"最高温度: {temp_max:.2f}°c")
print(f"中位数温度: {temp_median:.2f}°c")

# 计算分位数
q25 = np.percentile(temperature_data, 25)
q75 = np.percentile(temperature_data, 75)
print(f"第一四分位数: {q25:.2f}°c")
print(f"第三四分位数: {q75:.2f}°c")

# 识别极端天气事件
hot_days = temperature_data[temperature_data > temp_mean + 2 * temp_std]
cold_days = temperature_data[temperature_data < temp_mean - 2 * temp_std]

print(f"\n极端天气统计:")
print(f"异常热天数量: {len(hot_days)} 天")
print(f"异常冷天数量: {len(cold_days)} 天")

def analyze_data_distribution(data, title="数据分布分析"):
    """
    分析数据分布特性
    
    参数:
    data: 待分析的数据
    title: 分析标题
    """
    # 计算偏度和峰度
    from scipy import stats
    
    skewness = stats.skew(data)
    kurtosis = stats.kurtosis(data)
    
    print(f"{title}:")
    print(f"偏度: {skewness:.3f}")
    if skewness > 0:
        print("  - 数据右偏（正偏）")
    elif skewness < 0:
        print("  - 数据左偏（负偏）")
    else:
        print("  - 数据对称分布")
    
    print(f"峰度: {kurtosis:.3f}")
    if kurtosis > 0:
        print("  - 尖峰分布（比正态分布更尖锐）")
    elif kurtosis < 0:
        print("  - 平峰分布（比正态分布更平坦）")
    else:
        print("  - 正态峰度")
    
    # 进行正态性检验
    shapiro_stat, shapiro_p = stats.shapiro(data[:5000])  # shapiro-wilk检验（限制样本大小）
    print(f"shapiro-wilk正态性检验 p值: {shapiro_p:.6f}")
    if shapiro_p > 0.05:
        print("  - 数据符合正态分布（p > 0.05）")
    else:
        print("  - 数据不符合正态分布（p ≤ 0.05）")
    
    # 计算置信区间
    confidence_level = 0.95
    degrees_freedom = len(data) - 1
    sample_mean = np.mean(data)
    sample_standard_error = stats.sem(data)
    confidence_interval = stats.t.interval(confidence_level, degrees_freedom, 
                                         sample_mean, sample_standard_error)
    
    print(f"{confidence_level*100}% 置信区间: [{confidence_interval[0]:.3f}, {confidence_interval[1]:.3f}]")

# 分析不同类型的分布
print("=== 均匀分布分析 ===")
uniform_sample = np.random.uniform(-10, 10, 10000)
analyze_data_distribution(uniform_sample, "均匀分布")

print("\n=== 正态分布分析 ===")
normal_sample = np.random.normal(0, 1, 10000)
analyze_data_distribution(normal_sample, "正态分布")

print("\n=== 指数分布分析 ===")
exponential_sample = np.random.exponential(2, 10000)
analyze_data_distribution(exponential_sample, "指数分布")

def generate_correlated_data(n_samples=1000):
    """
    生成具有特定相关性的多维数据
    
    参数:
    n_samples: 样本数量
    
    返回:
    相关数据矩阵
    """
    # 定义相关系数矩阵
    correlation_matrix = np.array([
        [1.0, 0.8, 0.3],
        [0.8, 1.0, 0.1],
        [0.3, 0.1, 1.0]
    ])
    
    # 生成相关随机数
    mean = [0, 0, 0]
    data = np.random.multivariate_normal(mean, correlation_matrix, n_samples)
    
    return data

# 生成相关数据
correlated_data = generate_correlated_data(5000)

# 提取各列数据
var1 = correlated_data[:, 0]
var2 = correlated_data[:, 1]
var3 = correlated_data[:, 2]

# 计算皮尔逊相关系数
corr_12 = np.corrcoef(var1, var2)[0, 1]
corr_13 = np.corrcoef(var1, var3)[0, 1]
corr_23 = np.corrcoef(var2, var3)[0, 1]

print("变量间相关性分析:")
print(f"变量1与变量2的相关系数: {corr_12:.3f}")
print(f"变量1与变量3的相关系数: {corr_13:.3f}")
print(f"变量2与变量3的相关系数: {corr_23:.3f}")

# 计算协方差矩阵
covariance_matrix = np.cov(correlated_data.t)
print(f"\n协方差矩阵:")
print(covariance_matrix)

# 进行线性回归分析
from scipy import stats

slope, intercept, r_value, p_value, std_err = stats.linregress(var1, var2)
print(f"\n变量1对变量2的线性回归:")
print(f"斜率: {slope:.3f}")
print(f"截距: {intercept:.3f}")
print(f"相关系数r: {r_value:.3f}")
print(f"决定系数r²: {r_value**2:.3f}")
print(f"p值: {p_value:.6f}")

def perform_hypothesis_tests():
    """
    执行各种假设检验
    """
    # 生成两组样本数据
    sample1 = np.random.normal(100, 15, 100)  # 均值100，标准差15
    sample2 = np.random.normal(105, 15, 100)  # 均值105，标准差15
    
    print("=== 单样本t检验 ===")
    # 单样本t检验：检验样本均值是否等于特定值
    pop_mean = 100
    t_stat, p_val = stats.ttest_1samp(sample1, pop_mean)
    print(f"检验统计量: {t_stat:.3f}")
    print(f"p值: {p_val:.3f}")
    if p_val < 0.05:
        print("拒绝原假设：样本均值显著不同于100")
    else:
        print("接受原假设：样本均值与100无显著差异")
    
    print("\n=== 双样本t检验 ===")
    # 双样本t检验：检验两组样本的均值是否有显著差异
    t_stat, p_val = stats.ttest_ind(sample1, sample2)
    print(f"检验统计量: {t_stat:.3f}")
    print(f"p值: {p_val:.3f}")
    if p_val < 0.05:
        print("拒绝原假设：两组样本均值存在显著差异")
    else:
        print("接受原假设：两组样本均值无显著差异")
    
    print("\n=== 配对样本t检验 ===")
    # 配对样本t检验：检验配对观测值的差异
    before_treatment = np.random.normal(80, 10, 50)
    after_treatment = before_treatment + np.random.normal(5, 3, 50)  # 治疗效果提升5个单位
    t_stat, p_val = stats.ttest_rel(before_treatment, after_treatment)
    print(f"检验统计量: {t_stat:.3f}")
    print(f"p值: {p_val:.3f}")
    if p_val < 0.05:
        print("拒绝原假设：治疗前后存在显著差异")
    else:
        print("接受原假设：治疗前后无显著差异")
    
    print("\n=== 方差齐性检验 ===")
    # levene检验：检验方差齐性
    w_stat, p_val = stats.levene(sample1, sample2)
    print(f"levene检验统计量: {w_stat:.3f}")
    print(f"p值: {p_val:.3f}")
    if p_val < 0.05:
        print("拒绝原假设：两组样本方差不相等")
    else:
        print("接受原假设：两组样本方差相等")

perform_hypothesis_tests()

def simulate_survey(population_size=1000000, sample_size=1000, true_proportion=0.55):
    """
    模拟调查
    参数:
    population_size: 总体规模
    sample_size: 样本规模
    true_proportion: 真实支持率
    返回:
    调查结果
    """
    # 生成总体数据
    population = np.random.binomial(1, true_proportion, population_size)
    # 随机抽样
    sample_indices = np.random.choice(population_size, sample_size, replace=false)
    sample = population[sample_indices]
    # 计算样本比例
    sample_proportion = np.mean(sample)
    # 计算置信区间（使用正态近似）
    standard_error = np.sqrt(sample_proportion * (1 - sample_proportion) / sample_size)
    margin_of_error = 1.96 * standard_error  # 95%置信水平
    lower_bound = sample_proportion - margin_of_error
    upper_bound = sample_proportion + margin_of_error
    return {
        'sample_proportion': sample_proportion,
        'lower_bound': lower_bound,
        'upper_bound': upper_bound,
        'margin_of_error': margin_of_error,
        'true_proportion': true_proportion
    }
# 执行多次调查模拟
n_simulations = 100
results = []
for i in range(n_simulations):
    survey_result = simulate_survey()
    results.append(survey_result)
    # 检查真实值是否在置信区间内
    contains_true = (survey_result['lower_bound'] <= survey_result['true_proportion'] <= 
                     survey_result['upper_bound'])
    if i < 5:  # 只打印前5次结果
        print(f"调查 {i+1}: 支持率 {survey_result['sample_proportion']:.3f} "
              f"(95% ci: [{survey_result['lower_bound']:.3f}, {survey_result['upper_bound']:.3f}]) "
              f"包含真实值: {contains_true}")
# 计算置信区间的覆盖率
coverage_rate = sum(1 for result in results 
                   if result['lower_bound'] <= result['true_proportion'] <= result['upper_bound']) / n_simulations
print(f"\n=== 置信区间性能评估 ===")
print(f"理论覆盖率: 95%")
print(f"实际覆盖率: {coverage_rate*100:.1f}%")
print(f"模拟次数: {n_simulations}")
# 分析误差分布
sample_proportions = [result['sample_proportion'] for result in results]
errors = [prop - 0.55 for prop in sample_proportions]
mean_error = np.mean(errors)
std_error_dist = np.std(errors)
print(f"\n误差分析:")
print(f"平均误差: {mean_error:.4f}")
print(f"误差标准差: {std_error_dist:.4f}")

def generate_time_series(length=1000, trend=0.01, seasonality_period=50, noise_level=0.5):
    """
    生成时间序列数据
    
    参数:
    length: 序列长度
    trend: 趋势项
    seasonality_period: 季节性周期
    noise_level: 噪声水平
    
    返回:
    时间序列数据
    """
    time = np.arange(length)
    
    # 趋势项
    trend_component = trend * time
    
    # 季节性项
    seasonal_component = np.sin(2 * np.pi * time / seasonality_period)
    
    # 噪声项
    noise = np.random.normal(0, noise_level, length)
    
    # 自相关项（ar(1)过程）
    ar_coefficient = 0.8
    ar_process = np.zeros(length)
    ar_process[0] = noise[0]
    for i in range(1, length):
        ar_process[i] = ar_coefficient * ar_process[i-1] + noise[i]
    
    # 合成时间序列
    ts = trend_component + seasonal_component + ar_process
    
    return ts

# 生成时间序列数据
ts_data = generate_time_series(1000)

# 基本统计描述
print("=== 时间序列基本统计 ===")
print(f"均值: {np.mean(ts_data):.3f}")
print(f"标准差: {np.std(ts_data):.3f}")
print(f"最小值: {np.min(ts_data):.3f}")
print(f"最大值: {np.max(ts_data):.3f}")

# 移动平均
def moving_average(data, window_size):
    """计算移动平均"""
    weights = np.ones(window_size) / window_size
    return np.convolve(data, weights, mode='valid')

ma_10 = moving_average(ts_data, 10)
ma_50 = moving_average(ts_data, 50)

print(f"\n移动平均:")
print(f"10期移动平均最新值: {ma_10[-1]:.3f}")
print(f"50期移动平均最新值: {ma_50[-1]:.3f}")

# 自相关分析
def autocorrelation(data, max_lag=50):
    """计算自相关函数"""
    n = len(data)
    mean = np.mean(data)
    var = np.var(data)
    
    autocorr = []
    for lag in range(max_lag + 1):
        if lag == 0:
            autocorr.append(1.0)
        else:
            numerator = np.sum((data[:-lag] - mean) * (data[lag:] - mean))
            denominator = (n - lag) * var
            autocorr.append(numerator / denominator)
    
    return np.array(autocorr)

# 计算自相关
ac_values = autocorrelation(ts_data, 100)

print(f"\n自相关分析:")
print(f"滞后1自相关: {ac_values[1]:.3f}")
print(f"滞后10自相关: {ac_values[10]:.3f}")
print(f"滞后50自相关: {ac_values[50]:.3f}")

# 简单预测
def simple_forecast(data, forecast_steps=10):
    """基于最近值的简单预测"""
    last_value = data[-1]
    trend = data[-1] - data[-2]  # 最近的趋势
    forecast = [last_value + i * trend for i in range(1, forecast_steps + 1)]
    return np.array(forecast)

forecast_values = simple_forecast(ts_data, 20)
print(f"\n未来20期预测值范围: [{forecast_values[0]:.3f}, {forecast_values[-1]:.3f}]")

def assess_data_quality(data, name="dataset"):
    """
    评估数据质量
    
    参数:
    data: 待评估的数据
    name: 数据集名称
    """
    print(f"=== {name} 数据质量评估 ===")
    
    # 基本信息
    total_count = len(data)
    print(f"总记录数: {total_count}")
    
    # 缺失值检查
    missing_count = np.sum(np.isnan(data)) if data.dtype.kind in ['f', 'c'] else 0
    missing_percentage = missing_count / total_count * 100
    print(f"缺失值数量: {missing_count} ({missing_percentage:.2f}%)")
    
    # 异常值检测（使用iqr方法）
    q1 = np.percentile(data, 25)
    q3 = np.percentile(data, 75)
    iqr = q3 - q1
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr
    
    outliers = data[(data < lower_bound) | (data > upper_bound)]
    outlier_count = len(outliers)
    outlier_percentage = outlier_count / total_count * 100
    
    print(f"异常值数量: {outlier_count} ({outlier_percentage:.2f}%)")
    if outlier_count > 0:
        print(f"异常值范围: [{np.min(outliers):.3f}, {np.max(outliers):.3f}]")
    
    # 重复值检查
    unique_count = len(np.unique(data))
    duplicate_count = total_count - unique_count
    duplicate_percentage = duplicate_count / total_count * 100
    print(f"重复值数量: {duplicate_count} ({duplicate_percentage:.2f}%)")
    
    # 数据分布检查
    data_range = np.max(data) - np.min(data)
    print(f"数据范围: {data_range:.3f}")
    
    # 偏度和峰度
    from scipy import stats
    skewness = stats.skew(data)
    kurtosis = stats.kurtosis(data)
    print(f"偏度: {skewness:.3f}")
    print(f"峰度: {kurtosis:.3f}")
    
    print("-" * 40)

# 生成测试数据集
np.random.seed(42)

# 正常数据
normal_data = np.random.normal(50, 10, 1000)
assess_data_quality(normal_data, "正常分布数据")

# 包含异常值的数据
contaminated_data = np.concatenate([
    np.random.normal(50, 10, 950),
    np.random.normal(100, 5, 50)  # 添加一些异常值
])
assess_data_quality(contaminated_data, "包含异常值的数据")

# 包含缺失值的数据
incomplete_data = np.random.normal(50, 10, 1000)
missing_indices = np.random.choice(1000, 50, replace=false)
incomplete_data[missing_indices] = np.nan
assess_data_quality(incomplete_data, "包含缺失值的数据")

def simulate_ab_test(control_conversion_rate=0.1, treatment_conversion_rate=0.12, 
                     sample_size_per_group=10000, alpha=0.05):
    """
    模拟a/b测试
    
    参数:
    control_conversion_rate: 对照组转化率
    treatment_conversion_rate: 实验组转化率
    sample_size_per_group: 每组样本量
    alpha: 显著性水平
    
    返回:
    测试结果
    """
    # 生成对照组数据
    control_conversions = np.random.binomial(1, control_conversion_rate, sample_size_per_group)
    control_cr = np.mean(control_conversions)
    
    # 生成实验组数据
    treatment_conversions = np.random.binomial(1, treatment_conversion_rate, sample_size_per_group)
    treatment_cr = np.mean(treatment_conversions)
    
    # 计算差异
    difference = treatment_cr - control_cr
    pooled_cr = (np.sum(control_conversions) + np.sum(treatment_conversions)) / (2 * sample_size_per_group)
    
    # z检验
    se_diff = np.sqrt(pooled_cr * (1 - pooled_cr) * (1/sample_size_per_group + 1/sample_size_per_group))
    z_score = difference / se_diff
    p_value = 2 * (1 - stats.norm.cdf(abs(z_score)))  # 双尾检验
    
    # 置信区间
    margin_of_error = stats.norm.ppf(1 - alpha/2) * se_diff
    ci_lower = difference - margin_of_error
    ci_upper = difference + margin_of_error
    
    # 统计功效计算（简化版）
    effect_size = abs(treatment_conversion_rate - control_conversion_rate)
    power = stats.norm.cdf(
        (abs(difference) - stats.norm.ppf(1 - alpha/2) * se_diff) / se_diff
    ) if se_diff > 0 else 0
    
    return {
        'control_conversion_rate': control_cr,
        'treatment_conversion_rate': treatment_cr,
        'difference': difference,
        'z_score': z_score,
        'p_value': p_value,
        'significant': p_value < alpha,
        'ci_lower': ci_lower,
        'ci_upper': ci_upper,
        'power': power
    }

# 执行a/b测试模拟
print("=== a/b测试模拟结果 ===")

# 场景1：无实际差异
print("场景1：对照组转化率10%，实验组转化率10%（无差异）")
result1 = simulate_ab_test(0.1, 0.1, 10000)
print(f"对照组转化率: {result1['control_conversion_rate']:.4f}")
print(f"实验组转化率: {result1['treatment_conversion_rate']:.4f}")
print(f"差异: {result1['difference']:.4f}")
print(f"z值: {result1['z_score']:.3f}")
print(f"p值: {result1['p_value']:.6f}")
print(f"统计显著: {result1['significant']}")
print(f"95%置信区间: [{result1['ci_lower']:.4f}, {result1['ci_upper']:.4f}]")

print("\n场景2：存在实际差异")
result2 = simulate_ab_test(0.1, 0.12, 10000)
print(f"对照组转化率: {result2['control_conversion_rate']:.4f}")
print(f"实验组转化率: {result2['treatment_conversion_rate']:.4f}")
print(f"差异: {result2['difference']:.4f}")
print(f"z值: {result2['z_score']:.3f}")
print(f"p值: {result2['p_value']:.6f}")
print(f"统计显著: {result2['significant']}")
print(f"95%置信区间: [{result2['ci_lower']:.4f}, {result2['ci_upper']:.4f}]")

# 多次模拟以评估测试性能
def evaluate_ab_test_performance(true_difference, sample_size, n_simulations=1000):
    """评估a/b测试性能"""
    significant_results = 0
    
    for _ in range(n_simulations):
        result = simulate_ab_test(0.1, 0.1 + true_difference, sample_size)
        if result['significant']:
            significant_results += 1
    
    detection_rate = significant_results / n_simulations
    return detection_rate

print(f"\n=== a/b测试性能评估 ===")
print("检测真实差异的能力（统计功效）:")

differences = [0.005, 0.01, 0.02, 0.03]
sample_sizes = [5000, 10000, 20000]

for diff in differences:
    print(f"\n真实差异: {diff}")
    for size in sample_sizes:
        power = evaluate_ab_test_performance(diff, size, 500)
        print(f"  样本量 {size}: 检测率 {power:.1%}")

def simple_linear_regression(x, y):
    """
    简单线性回归
    
    参数:
    x: 自变量
    y: 因变量
    
    返回:
    回归参数和统计信息
    """
    # 计算回归系数
    n = len(x)
    x_mean = np.mean(x)
    y_mean = np.mean(y)
    
    # 斜率和截距
    slope = np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean)**2)
    intercept = y_mean - slope * x_mean
    
    # 预测值
    y_pred = slope * x + intercept
    
    # 残差
    residuals = y - y_pred
    
    # 统计量
    ss_res = np.sum(residuals**2)
    ss_tot = np.sum((y - y_mean)**2)
    r_squared = 1 - (ss_res / ss_tot)
    
    # 标准误差
    mse = ss_res / (n - 2)
    se_slope = np.sqrt(mse / np.sum((x - x_mean)**2))
    se_intercept = se_slope * np.sqrt(np.sum(x**2) / n)
    
    # t统计量和p值
    t_slope = slope / se_slope
    t_intercept = intercept / se_intercept
    
    # 简化的p值计算（假设大样本）
    from scipy import stats
    p_slope = 2 * (1 - stats.t.cdf(abs(t_slope), n - 2))
    p_intercept = 2 * (1 - stats.t.cdf(abs(t_intercept), n - 2))
    
    return {
        'slope': slope,
        'intercept': intercept,
        'r_squared': r_squared,
        'se_slope': se_slope,
        'se_intercept': se_intercept,
        't_slope': t_slope,
        't_intercept': t_intercept,
        'p_slope': p_slope,
        'p_intercept': p_intercept,
        'predictions': y_pred,
        'residuals': residuals
    }

# 生成带噪声的线性关系数据
np.random.seed(42)
x_data = np.linspace(0, 10, 100)
true_slope = 2.5
true_intercept = 1.0
noise = np.random.normal(0, 2, 100)
y_data = true_slope * x_data + true_intercept + noise

# 执行回归分析
regression_result = simple_linear_regression(x_data, y_data)

print("=== 线性回归分析结果 ===")
print(f"估计斜率: {regression_result['slope']:.3f} (真实值: {true_slope})")
print(f"估计截距: {regression_result['intercept']:.3f} (真实值: {true_intercept})")
print(f"决定系数 r²: {regression_result['r_squared']:.3f}")
print(f"斜率标准误差: {regression_result['se_slope']:.3f}")
print(f"截距标准误差: {regression_result['se_intercept']:.3f}")
print(f"斜率t统计量: {regression_result['t_slope']:.3f}")
print(f"斜率p值: {regression_result['p_slope']:.6f}")
print(f"截距p值: {regression_result['p_intercept']:.6f}")

# 模型诊断
residuals = regression_result['residuals']
print(f"\n残差分析:")
print(f"残差均值: {np.mean(residuals):.6f} (应接近0)")
print(f"残差标准差: {np.std(residuals):.3f}")

# 检查残差的正态性
from scipy import stats
shapiro_stat, shapiro_p = stats.shapiro(residuals[:5000])  # 限制样本大小
print(f"残差正态性检验 p值: {shapiro_p:.6f}")

def monte_carlo_pi_estimation(n_samples=1000000):
    """
    使用蒙特卡洛方法估算π值
    
    参数:
    n_samples: 样本数量
    
    返回:
    π的估算值
    """
    # 在单位正方形内生成随机点
    x = np.random.uniform(-1, 1, n_samples)
    y = np.random.uniform(-1, 1, n_samples)
    
    # 计算点到原点的距离
    distances = np.sqrt(x**2 + y**2)
    
    # 统计在单位圆内的点数
    points_in_circle = np.sum(distances <= 1)
    
    # 估算π值
    pi_estimate = 4 * points_in_circle / n_samples
    
    return pi_estimate

def monte_carlo_integration(func, a, b, n_samples=1000000):
    """
    使用蒙特卡洛方法进行数值积分
    
    参数:
    func: 被积函数
    a, b: 积分区间
    n_samples: 样本数量
    
    返回:
    积分估算值
    """
    # 在区间[a,b]内生成随机点
    x = np.random.uniform(a, b, n_samples)
    
    # 计算函数值
    y = func(x)
    
    # 估算积分值
    integral_estimate = (b - a) * np.mean(y)
    
    return integral_estimate

# 估算π值
pi_estimate = monte_carlo_pi_estimation(1000000)
print(f"=== 蒙特卡洛π值估算 ===")
print(f"估算值: {pi_estimate:.6f}")
print(f"真实值: {np.pi:.6f}")
print(f"误差: {abs(pi_estimate - np.pi):.6f}")

# 数值积分示例：计算sin(x)在[0,π]上的积分
def test_function(x):
    return np.sin(x)

integration_result = monte_carlo_integration(test_function, 0, np.pi, 1000000)
analytical_result = 2  # sin(x)在[0,π]上的定积分

print(f"\n=== 蒙特卡洛数值积分 ===")
print(f"sin(x)在[0,π]上的积分:")
print(f"蒙特卡洛估算: {integration_result:.6f}")
print(f"解析解: {analytical_result:.6f}")
print(f"误差: {abs(integration_result - analytical_result):.6f}")

# 风险评估示例
def monte_carlo_risk_assessment(initial_investment=10000, years=10, n_simulations=10000):
    """
    使用蒙特卡洛方法进行投资风险评估
    
    参数:
    initial_investment: 初始投资金额
    years: 投资年限
    n_simulations: 模拟次数
    
    返回:
    投资结果统计
    """
    # 假设年收益率服从正态分布（均值8%，标准差15%）
    annual_returns = np.random.normal(0.08, 0.15, (years, n_simulations))
    
    # 计算累计收益
    cumulative_returns = np.ones((years + 1, n_simulations))
    for year in range(1, years + 1):
        cumulative_returns[year] = cumulative_returns[year-1] * (1 + annual_returns[year-1])
    
    # 最终价值
    final_values = initial_investment * cumulative_returns[-1]
    
    # 统计分析
    mean_final_value = np.mean(final_values)
    median_final_value = np.median(final_values)
    std_final_value = np.std(final_values)
    
    # 计算损失概率（最终价值低于初始投资）
    loss_probability = np.sum(final_values < initial_investment) / n_simulations
    
    # 计算var（value at risk）
    var_95 = np.percentile(final_values, 5)
    var_99 = np.percentile(final_values, 1)
    
    return {
        'mean_final_value': mean_final_value,
        'median_final_value': median_final_value,
        'std_final_value': std_final_value,
        'loss_probability': loss_probability,
        'var_95': var_95,
        'var_99': var_99,
        'final_values': final_values
    }

# 执行风险评估
risk_result = monte_carlo_risk_assessment(10000, 10, 10000)

print(f"\n=== 投资风险评估结果 ===")
print(f"初始投资: $10,000")
print(f"投资期限: 10年")
print(f"平均最终价值: ${risk_result['mean_final_value']:,.2f}")
print(f"中位数最终价值: ${risk_result['median_final_value']:,.2f}")
print(f"标准差: ${risk_result['std_final_value']:,.2f}")
print(f"亏损概率: {risk_result['loss_probability']:.1%}")
print(f"95% var: ${risk_result['var_95']:,.2f}")
print(f"99% var: ${risk_result['var_99']:,.2f}")