Python实例题：Python计算数理统计

于 2025-06-08 19:00:00 发布
阅读量489 收藏 6
点赞数 3
CC 4.0 BY-SA版权
分类专栏： 实例 文章标签： python 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/shangzhiqi/article/details/148502732

目录

Python实例题

题目

代码实现

实现原理

描述性统计：

参数估计：

假设检验：

相关与回归分析：

可视化功能：

关键代码解析

1. 描述性统计

2. 置信区间计算

3. 假设检验

4. 线性回归分析

使用说明

安装依赖：

基本用法：

示例输出：

扩展建议

增强功能：

用户界面：

性能优化：

教学辅助：

Python实例题

题目

Python计算数理统计

代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
import pandas as pd
from sympy import symbols, simplify, latex
from statsmodels.formula.api import ols
import statsmodels.api as sm

class StatisticalAnalysis:
    """数理统计分析类，支持描述性统计、参数估计和假设检验"""
    
    def __init__(self):
        """初始化分析器"""
        pass
    
    def descriptive_statistics(self, data):
        """
        计算描述性统计量
        
        参数:
            data: 样本数据（列表或NumPy数组）
            
        返回:
            dict: 包含各种描述性统计量的字典
        """
        results = {
            '样本大小': len(data),
            '均值': np.mean(data),
            '中位数': np.median(data),
            '众数': stats.mode(data, keepdims=True)[0][0],
            '标准差': np.std(data, ddof=1),  # 无偏估计
            '方差': np.var(data, ddof=1),    # 无偏估计
            '最小值': np.min(data),
            '最大值': np.max(data),
            '极差': np.max(data) - np.min(data),
            '四分位数Q1': np.percentile(data, 25),
            '四分位数Q2': np.percentile(data, 50),
            '四分位数Q3': np.percentile(data, 75),
            '四分位距': np.percentile(data, 75) - np.percentile(data, 25),
            '偏度': stats.skew(data),
            '峰度': stats.kurtosis(data, fisher=True)  # Fisher校正的峰度
        }
        return results
    
    def confidence_interval(self, data, confidence_level=0.95):
        """
        计算均值的置信区间
        
        参数:
            data: 样本数据（列表或NumPy数组）
            confidence_level: 置信水平，默认为0.95
            
        返回:
            tuple: (置信下限, 置信上限)
        """
        n = len(data)
        mean = np.mean(data)
        std = np.std(data, ddof=1)  # 样本标准差
        alpha = 1 - confidence_level
        
        # 计算t分布的临界值
        t_critical = stats.t.ppf(1 - alpha/2, df=n-1)
        
        # 计算标准误差
        standard_error = std / np.sqrt(n)
        
        # 计算置信区间
        lower = mean - t_critical * standard_error
        upper = mean + t_critical * standard_error
        
        return (lower, upper)
    
    def hypothesis_test_one_sample(self, data, mu0, alternative='two-sided', alpha=0.05):
        """
        单样本t检验
        
        参数:
            data: 样本数据（列表或NumPy数组）
            mu0: 原假设的总体均值
            alternative: 备择假设方向，'two-sided'（双侧）、'greater'（大于）或'less'（小于）
            alpha: 显著性水平，默认为0.05
            
        返回:
            dict: 包含检验统计量、p值和检验结论的字典
        """
        # 执行t检验
        t_stat, p_value = stats.ttest_1samp(data, mu0)
        
        # 根据备择假设调整p值
        if alternative == 'greater' and t_stat > 0:
            p_value = p_value / 2
        elif alternative == 'less' and t_stat < 0:
            p_value = p_value / 2
        elif alternative in ['greater', 'less']:
            p_value = 1 - p_value / 2
        
        # 得出结论
        reject_null = p_value < alpha
        conclusion = f"在显著性水平 {alpha} 下，{'拒绝' if reject_null else '不拒绝'} 原假设"
        
        results = {
            '检验统计量': t_stat,
            'p值': p_value,
            '自由度': len(data) - 1,
            '结论': conclusion
        }
        
        return results
    
    def hypothesis_test_two_samples(self, data1, data2, equal_var=True, alternative='two-sided', alpha=0.05):
        """
        两样本t检验
        
        参数:
            data1: 第一个样本数据
            data2: 第二个样本数据
            equal_var: 是否假设两总体方差相等，默认为True
            alternative: 备择假设方向，'two-sided'（双侧）、'greater'（大于）或'less'（小于）
            alpha: 显著性水平，默认为0.05
            
        返回:
            dict: 包含检验统计量、p值和检验结论的字典
        """
        # 执行t检验
        t_stat, p_value = stats.ttest_ind(data1, data2, equal_var=equal_var)
        
        # 根据备择假设调整p值
        if alternative == 'greater' and t_stat > 0:
            p_value = p_value / 2
        elif alternative == 'less' and t_stat < 0:
            p_value = p_value / 2
        elif alternative in ['greater', 'less']:
            p_value = 1 - p_value / 2
        
        # 得出结论
        reject_null = p_value < alpha
        conclusion = f"在显著性水平 {alpha} 下，{'拒绝' if reject_null else '不拒绝'} 原假设"
        
        results = {
            '检验统计量': t_stat,
            'p值': p_value,
            '自由度': len(data1) + len(data2) - 2 if equal_var else None,
            '结论': conclusion
        }
        
        return results
    
    def anova_one_way(self, *groups):
        """
        单因素方差分析
        
        参数:
            *groups: 多个样本组（每个组是一个列表或NumPy数组）
            
        返回:
            dict: 包含F统计量、p值和检验结论的字典
        """
        # 执行单因素方差分析
        f_stat, p_value = stats.f_oneway(*groups)
        
        # 得出结论
        alpha = 0.05
        reject_null = p_value < alpha
        conclusion = f"在显著性水平 {alpha} 下，{'至少有一组均值存在显著差异' if reject_null else '各组均值无显著差异'}"
        
        results = {
            'F统计量': f_stat,
            'p值': p_value,
            '结论': conclusion
        }
        
        return results
    
    def correlation_analysis(self, X, Y, method='pearson'):
        """
        相关性分析
        
        参数:
            X: 第一个变量的数据
            Y: 第二个变量的数据
            method: 相关系数计算方法，'pearson'（皮尔逊）或'spearman'（斯皮尔曼）
            
        返回:
            dict: 包含相关系数、p值和结论的字典
        """
        if method == 'pearson':
            corr, p_value = stats.pearsonr(X, Y)
        elif method == 'spearman':
            corr, p_value = stats.spearmanr(X, Y)
        else:
            raise ValueError("不支持的相关系数计算方法")
        
        # 得出结论
        alpha = 0.05
        significant = p_value < alpha
        strength = ''
        if abs(corr) < 0.2:
            strength = '极弱相关'
        elif abs(corr) < 0.4:
            strength = '弱相关'
        elif abs(corr) < 0.6:
            strength = '中等相关'
        elif abs(corr) < 0.8:
            strength = '强相关'
        else:
            strength = '极强相关'
        
        direction = '正相关' if corr > 0 else '负相关'
        
        results = {
            f'{method}相关系数': corr,
            'p值': p_value,
            '显著性': f"在显著性水平 {alpha} 下，{'显著' if significant else '不显著'}",
            '相关强度和方向': f"{strength}，{direction}"
        }
        
        return results
    
    def linear_regression(self, X, Y):
        """
        简单线性回归分析
        
        参数:
            X: 自变量数据
            Y: 因变量数据
            
        返回:
            dict: 包含回归系数、截距、R²值和p值的字典
        """
        # 添加常数项
        X_with_const = sm.add_constant(X)
        
        # 执行线性回归
        model = sm.OLS(Y, X_with_const).fit()
        
        # 提取结果
        results = {
            '截距': model.params[0],
            '斜率': model.params[1],
            'R²值': model.rsquared,
            '调整R²值': model.rsquared_adj,
            'F统计量': model.fvalue,
            'F检验p值': model.f_pvalue,
            '斜率p值': model.pvalues[1],
            '回归方程': f"Y = {model.params[0]:.4f} + {model.params[1]:.4f}X",
            '模型摘要': model.summary()
        }
        
        return results
    
    def chi_square_test(self, observed):
        """
        卡方检验（拟合优度检验或独立性检验）
        
        参数:
            observed: 观察频数的二维数组或列表的列表
            
        返回:
            dict: 包含卡方统计量、p值、自由度和结论的字典
        """
        # 执行卡方检验
        chi2, p, dof, expected = stats.chi2_contingency(observed)
        
        # 得出结论
        alpha = 0.05
        reject_null = p < alpha
        conclusion = f"在显著性水平 {alpha} 下，{'拒绝' if reject_null else '不拒绝'} 原假设"
        
        results = {
            '卡方统计量': chi2,
            'p值': p,
            '自由度': dof,
            '期望频数': expected,
            '结论': conclusion
        }
        
        return results
    
    def plot_histogram(self, data, bins=30, title=None):
        """
        绘制样本直方图
        
        参数:
            data: 样本数据
            bins: 直方图的箱数，默认为30
            title: 图像标题，默认为None
        """
        plt.figure(figsize=(10, 6))
        sns.histplot(data, bins=bins, kde=True)
        
        plt.grid(True)
        plt.axvline(x=np.mean(data), color='r', linestyle='--', label=f'均值: {np.mean(data):.2f}')
        plt.axvline(x=np.median(data), color='g', linestyle='--', label=f'中位数: {np.median(data):.2f}')
        
        if title:
            plt.title(title, fontsize=14)
        else:
            plt.title(f"样本直方图", fontsize=14)
        
        plt.xlabel('值', fontsize=12)
        plt.ylabel('频率', fontsize=12)
        plt.legend()
        plt.show()
    
    def plot_boxplot(self, data, title=None):
        """
        绘制箱线图
        
        参数:
            data: 样本数据（可以是单个列表或多个列表的列表）
            title: 图像标题，默认为None
        """
        plt.figure(figsize=(10, 6))
        
        if isinstance(data[0], list) or isinstance(data[0], np.ndarray):
            # 多个数据集
            plt.boxplot(data)
            plt.xticks(range(1, len(data) + 1), [f"数据{i+1}" for i in range(len(data))])
        else:
            # 单个数据集
            plt.boxplot(data)
        
        plt.grid(True)
        
        if title:
            plt.title(title, fontsize=14)
        else:
            plt.title(f"箱线图", fontsize=14)
        
        plt.ylabel('值', fontsize=12)
        plt.show()
    
    def plot_scatter(self, X, Y, title=None, xlabel=None, ylabel=None):
        """
        绘制散点图
        
        参数:
            X: 自变量数据
            Y: 因变量数据
            title: 图像标题，默认为None
            xlabel: x轴标签，默认为None
            ylabel: y轴标签，默认为None
        """
        plt.figure(figsize=(10, 6))
        plt.scatter(X, Y, alpha=0.7)
        
        # 计算并绘制回归线
        slope, intercept, r_value, p_value, std_err = stats.linregress(X, Y)
        plt.plot(X, intercept + slope*X, 'r-', label=f'回归线: y = {intercept:.4f} + {slope:.4f}x')
        
        plt.grid(True)
        
        if title:
            plt.title(title, fontsize=14)
        else:
            plt.title(f"散点图 (r={r_value:.4f})", fontsize=14)
        
        if xlabel:
            plt.xlabel(xlabel, fontsize=12)
        else:
            plt.xlabel('X', fontsize=12)
        
        if ylabel:
            plt.ylabel(ylabel, fontsize=12)
        else:
            plt.ylabel('Y', fontsize=12)
        
        plt.legend()
        plt.show()


# 示例使用
def example_usage():
    analysis = StatisticalAnalysis()
    
    print("\n===== 描述性统计示例 =====")
    # 生成样本数据
    np.random.seed(42)  # 设置随机种子，确保结果可重现
    data = np.random.normal(loc=50, scale=10, size=100)
    
    # 计算描述性统计量
    stats = analysis.descriptive_statistics(data)
    for key, value in stats.items():
        print(f"{key}: {value:.4f}" if isinstance(value, float) else f"{key}: {value}")
    
    # 绘制直方图
    analysis.plot_histogram(data, title="样本数据直方图")
    # 绘制箱线图
    analysis.plot_boxplot(data, title="样本数据箱线图")
    
    print("\n===== 参数估计示例 =====")
    # 计算置信区间
    lower, upper = analysis.confidence_interval(data, confidence_level=0.95)
    print(f"95%置信区间: [{lower:.4f}, {upper:.4f}]")
    
    print("\n===== 假设检验示例 =====")
    # 单样本t检验
    test_results = analysis.hypothesis_test_one_sample(data, mu0=50)
    for key, value in test_results.items():
        print(f"{key}: {value}")
    
    # 两样本t检验
    data2 = np.random.normal(loc=55, scale=10, size=100)
    test_results = analysis.hypothesis_test_two_samples(data, data2, equal_var=True)
    for key, value in test_results.items():
        print(f"{key}: {value}")
    
    # 单因素方差分析
    group1 = np.random.normal(loc=50, scale=10, size=50)
    group2 = np.random.normal(loc=55, scale=10, size=50)
    group3 = np.random.normal(loc=60, scale=10, size=50)
    anova_results = analysis.anova_one_way(group1, group2, group3)
    for key, value in anova_results.items():
        print(f"{key}: {value}")
    
    print("\n===== 相关性分析示例 =====")
    # 生成相关数据
    X = np.random.normal(0, 1, 100)
    Y = 2*X + np.random.normal(0, 1, 100)
    
    # 计算相关系数
    corr_results = analysis.correlation_analysis(X, Y, method='pearson')
    for key, value in corr_results.items():
        print(f"{key}: {value}")
    
    # 绘制散点图
    analysis.plot_scatter(X, Y, title="X和Y的散点图", xlabel="X", ylabel="Y")
    
    print("\n===== 线性回归分析示例 =====")
    # 执行线性回归
    regression_results = analysis.linear_regression(X, Y)
    for key, value in regression_results.items():
        if key != '模型摘要':
            print(f"{key}: {value}")
    
    print("\n===== 卡方检验示例 =====")
    # 创建列联表
    observed = np.array([[20, 30, 10], [15, 25, 20]])
    chi2_results = analysis.chi_square_test(observed)
    for key, value in chi2_results.items():
        if key != '期望频数':
            print(f"{key}: {value}")
        else:
            print(f"{key}:\n{value}")


if __name__ == "__main__":
    example_usage()

实现原理

这个数理统计计算工具基于以下技术实现：

  • 描述性统计：

    • 计算均值、中位数、众数、标准差等基本统计量
    • 分析数据分布特征（偏度、峰度）
    • 提供直观的数据可视化（直方图、箱线图）

  • 参数估计：

    • 计算样本均值的置信区间
    • 基于 t 分布进行区间估计
    • 支持不同置信水平的选择

  • 假设检验：

    • 实现单样本和两样本 t 检验
    • 支持单因素方差分析（ANOVA）
    • 提供卡方检验（拟合优度和独立性检验）

  • 相关与回归分析：

    • 计算皮尔逊和斯皮尔曼相关系数
    • 执行简单线性回归分析
    • 评估回归模型的拟合优度和显著性

  • 可视化功能：

    • 使用 Matplotlib 和 Seaborn 绘制统计图表
    • 直观展示数据分布和变量关系
    • 支持自定义图表样式和标题

关键代码解析

1. 描述性统计

def descriptive_statistics(self, data):
    """计算描述性统计量"""
    results = {
        '样本大小': len(data),
        '均值': np.mean(data),
        '中位数': np.median(data),
        '众数': stats.mode(data, keepdims=True)[0][0],
        '标准差': np.std(data, ddof=1),
        '方差': np.var(data, ddof=1),
        '最小值': np.min(data),
        '最大值': np.max(data),
        '极差': np.max(data) - np.min(data),
        '四分位数Q1': np.percentile(data, 25),
        '四分位数Q2': np.percentile(data, 50),
        '四分位数Q3': np.percentile(data, 75),
        '四分位距': np.percentile(data, 75) - np.percentile(data, 25),
        '偏度': stats.skew(data),
        '峰度': stats.kurtosis(data, fisher=True)
    }
    return results

2. 置信区间计算

def confidence_interval(self, data, confidence_level=0.95):
    """计算均值的置信区间"""
    n = len(data)
    mean = np.mean(data)
    std = np.std(data, ddof=1)
    alpha = 1 - confidence_level
    
    t_critical = stats.t.ppf(1 - alpha/2, df=n-1)
    standard_error = std / np.sqrt(n)
    
    lower = mean - t_critical * standard_error
    upper = mean + t_critical * standard_error
    
    return (lower, upper)

3. 假设检验

def hypothesis_test_one_sample(self, data, mu0, alternative='two-sided', alpha=0.05):
    """单样本t检验"""
    t_stat, p_value = stats.ttest_1samp(data, mu0)
    
    if alternative == 'greater' and t_stat > 0:
        p_value = p_value / 2
    elif alternative == 'less' and t_stat < 0:
        p_value = p_value / 2
    elif alternative in ['greater', 'less']:
        p_value = 1 - p_value / 2
    
    reject_null = p_value < alpha
    conclusion = f"在显著性水平 {alpha} 下，{'拒绝' if reject_null else '不拒绝'} 原假设"
    
    results = {
        '检验统计量': t_stat,
        'p值': p_value,
        '自由度': len(data) - 1,
        '结论': conclusion
    }
    
    return results

4. 线性回归分析

def linear_regression(self, X, Y):
    """简单线性回归分析"""
    X_with_const = sm.add_constant(X)
    model = sm.OLS(Y, X_with_const).fit()
    
    results = {
        '截距': model.params[0],
        '斜率': model.params[1],
        'R²值': model.rsquared,
        '调整R²值': model.rsquared_adj,
        'F统计量': model.fvalue,
        'F检验p值': model.f_pvalue,
        '斜率p值': model.pvalues[1],
        '回归方程': f"Y = {model.params[0]:.4f} + {model.params[1]:.4f}X",
        '模型摘要': model.summary()
    }
    
    return results

使用说明

  • 安装依赖：

pip install numpy matplotlib scipy seaborn pandas statsmodels

  • 基本用法：

from statistical_analysis import StatisticalAnalysis

# 创建分析器实例
analysis = StatisticalAnalysis()

# 生成样本数据
data = np.random.normal(loc=50, scale=10, size=100)

# 计算描述性统计量
stats = analysis.descriptive_statistics(data)
for key, value in stats.items():
    print(f"{key}: {value}")

# 计算置信区间
lower, upper = analysis.confidence_interval(data, confidence_level=0.95)
print(f"95%置信区间: [{lower:.4f}, {upper:.4f}]")

# 单样本t检验
test_results = analysis.hypothesis_test_one_sample(data, mu0=50)
for key, value in test_results.items():
    print(f"{key}: {value}")

# 相关性分析
X = np.random.normal(0, 1, 100)
Y = 2*X + np.random.normal(0, 1, 100)
corr_results = analysis.correlation_analysis(X, Y)
for key, value in corr_results.items():
    print(f"{key}: {value}")

# 线性回归分析
regression_results = analysis.linear_regression(X, Y)
for key, value in regression_results.items():
    if key != '模型摘要':
        print(f"{key}: {value}")

# 绘制散点图
analysis.plot_scatter(X, Y, title="X和Y的散点图")

  • 示例输出：

样本大小: 100
均值: 49.8323
中位数: 49.7734
众数: 49.2435
标准差: 9.9247
方差: 98.4995
最小值: 27.2072
最大值: 74.4235
极差: 47.2163
四分位数Q1: 42.7174
四分位数Q2: 49.7734
四分位数Q3: 57.1504
四分位距: 14.4330
偏度: -0.0423
峰度: -0.2341
95%置信区间: [47.8743, 51.7903]
检验统计量: -0.1694
p值: 0.8652
自由度: 99
结论: 在显著性水平 0.05 下，不拒绝原假设
pearson相关系数: 0.8742
p值: 3.1643e-26
显著性: 在显著性水平 0.05 下，显著
相关强度和方向: 强相关，正相关
截距: -0.0616
斜率: 1.9774
R²值: 0.7642
调整R²值: 0.7619
F统计量: 331.1444
F检验p值: 3.16428e-26
斜率p值: 3.16428e-26
回归方程: Y = -0.0616 + 1.9774X

扩展建议

  • 增强功能：

    • 添加多元线性回归分析
    • 实现时间序列分析功能
    • 增加非参数检验方法
    • 支持贝叶斯统计分析

  • 用户界面：

    • 开发命令行交互界面
    • 创建图形界面（如使用 Tkinter 或 PyQt）
    • 实现 Web 界面（如使用 Flask 或 Django）

  • 性能优化：

    • 针对大规模数据计算进行优化
    • 添加并行计算支持
    • 优化内存使用

  • 教学辅助：

    • 添加统计概念解释和公式推导
    • 提供假设检验决策规则说明
    • 实现交互式可视化（如动态调整参数）
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值