Python实战用NumPy和SciPy玩转折叠正态分布附完整代码当你处理传感器读数、金融波动分析或信号处理时经常会遇到只关注数值大小而忽略正负的场景。这时候折叠正态分布Folded Normal Distribution就派上了大用场。它本质上是将普通正态分布沿y轴折叠只保留非负部分的数据分布。1. 理解折叠正态分布的核心概念折叠正态分布描述的是对正态随机变量取绝对值后的分布形态。假设原始变量Y服从N(μ,σ²)那么X|Y|就形成了折叠正态分布。这种分布在以下场景特别有用工业质检零件尺寸偏差的绝对值分析金融分析价格波动幅度的建模信号处理振幅强度的概率分布研究与普通正态分布相比折叠正态分布有几点关键差异特性正态分布折叠正态分布定义域(-∞, ∞)[0, ∞)对称性关于μ对称非对称μ≠0时峰值数量单峰可能双峰当μ远离0时注意当μ0时折叠正态分布退化为半正态分布这是工程中最常见的特例。2. 用NumPy生成折叠正态数据让我们从实战开始用NumPy生成折叠正态分布的样本数据import numpy as np import matplotlib.pyplot as plt def generate_folded_normal(mu, sigma, size1000): 生成折叠正态分布样本 参数 mu: 原始正态分布的均值 sigma: 原始正态分布的标准差 size: 样本数量 返回 折叠后的样本数组 normal_samples np.random.normal(mu, sigma, size) return np.abs(normal_samples) # 生成三种不同参数的折叠正态分布 samples1 generate_folded_normal(0, 1, 10000) # 半正态情况 samples2 generate_folded_normal(1, 1, 10000) # 中等偏移 samples3 generate_folded_normal(3, 1, 10000) # 显著偏移可视化这些分布能直观理解参数影响plt.figure(figsize(12, 6)) bins np.linspace(0, 6, 50) plt.hist(samples1, bins, alpha0.5, labelμ0, σ1) plt.hist(samples2, bins, alpha0.5, labelμ1, σ1) plt.hist(samples3, bins, alpha0.5, labelμ3, σ1) plt.title(不同参数下的折叠正态分布对比) plt.xlabel(数值) plt.ylabel(频数) plt.legend() plt.grid(True) plt.show()运行这段代码你会观察到当μ0时分布从0单调递减随着μ增大分布逐渐右移并可能出现双峰σ控制着分布的扩散程度3. SciPy中的精确概率计算虽然NumPy适合生成随机样本但SciPy提供了更精确的概率计算工具。我们可以定义折叠正态分布的PDF和CDFfrom scipy.stats import norm from scipy.integrate import quad def folded_normal_pdf(x, mu, sigma): 折叠正态分布的概率密度函数 if x 0: return 0 return (norm.pdf(x, mu, sigma) norm.pdf(-x, mu, sigma)) def folded_normal_cdf(x, mu, sigma): 折叠正态分布的累积分布函数 if x 0: return 0 return norm.cdf(x, mu, sigma) - norm.cdf(-x, mu, sigma)验证这些函数的正确性x_vals np.linspace(0, 5, 500) pdf_vals [folded_normal_pdf(x, 1, 1) for x in x_vals] cdf_vals [folded_normal_cdf(x, 1, 1) for x in x_vals] plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.plot(x_vals, pdf_vals, r-, lw2) plt.title(PDF (μ1, σ1)) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(x_vals, cdf_vals, b-, lw2) plt.title(CDF (μ1, σ1)) plt.grid(True) plt.tight_layout() plt.show()对于更复杂的统计量计算比如均值和方差def folded_normal_mean(mu, sigma): 计算折叠正态分布的期望 term1 mu * (1 - 2 * norm.cdf(-mu/sigma)) term2 sigma * np.sqrt(2/np.pi) * np.exp(-mu**2/(2*sigma**2)) return term1 term2 def folded_normal_var(mu, sigma): 计算折叠正态分布的方差 EX folded_normal_mean(mu, sigma) EX2 mu**2 sigma**2 return EX2 - EX**2 # 计算示例 print(fμ1, σ1时的期望: {folded_normal_mean(1, 1):.4f}) print(fμ1, σ1时的方差: {folded_normal_var(1, 1):.4f})4. 实际应用案例质量控制分析假设我们有一批零件的长度设计值为10mm允许有±0.2mm的误差。实际测量发现误差服从N(0.05, 0.1²)的正态分布。我们想分析绝对误差的分布情况。# 参数设置 design_length 10 # mm mu_error 0.05 # mm sigma_error 0.1 # mm # 生成模拟数据 errors np.random.normal(mu_error, sigma_error, 1000) abs_errors np.abs(errors) # 关键指标计算 prob_within_tolerance folded_normal_cdf(0.2, mu_error, sigma_error) expected_abs_error folded_normal_mean(mu_error, sigma_error) print(f绝对误差在允许范围内的概率: {prob_within_tolerance:.2%}) print(f预期绝对误差: {expected_abs_error:.4f} mm) # 可视化 plt.figure(figsize(10, 5)) x np.linspace(0, 0.3, 100) plt.plot(x, [folded_normal_pdf(xi, mu_error, sigma_error) for xi in x], r-, lw2, label理论PDF) plt.hist(abs_errors, bins30, densityTrue, alpha0.6, label模拟数据) plt.axvline(0.2, colork, linestyle--, label允许上限) plt.title(零件长度绝对误差分布) plt.xlabel(绝对误差(mm)) plt.ylabel(概率密度) plt.legend() plt.grid(True) plt.show()这个案例展示了如何将折叠正态分布应用于实际质量控制场景。通过计算我们能够量化产品合格概率和预期误差水平为生产决策提供数据支持。5. 高级技巧参数估计与模型拟合当我们有实际观测数据时如何估计折叠正态分布的参数这里介绍两种方法方法一矩估计法def estimate_params_moment(data): 通过矩估计法估计μ和σ m1 np.mean(data) m2 np.mean(data**2) # 解非线性方程组 from scipy.optimize import root def equations(p): mu, sigma p E_X folded_normal_mean(mu, sigma) E_X2 m2 # μ² σ² return [E_X - m1, mu**2 sigma**2 - E_X2] sol root(equations, [1, 1]) return sol.x[0], sol.x[1]方法二最大似然估计from scipy.optimize import minimize def folded_normal_loglik(params, data): 折叠正态分布的负对数似然函数 mu, sigma params if sigma 0: return np.inf pdf_vals [folded_normal_pdf(x, mu, sigma) for x in data] return -np.sum(np.log(pdf_vals)) def estimate_params_mle(data): 最大似然估计 # 使用矩估计结果作为初始值 mu_init, sigma_init estimate_params_moment(data) res minimize(folded_normal_loglik, [mu_init, sigma_init], args(data,), bounds[(None, None), (1e-6, None)]) return res.x[0], res.x[1]比较两种方法的估计效果# 生成测试数据 true_mu, true_sigma 1.5, 0.8 test_data generate_folded_normal(true_mu, true_sigma, 1000) # 参数估计 mu_moment, sigma_moment estimate_params_moment(test_data) mu_mle, sigma_mle estimate_params_mle(test_data) print(f真实参数: μ{true_mu:.2f}, σ{true_sigma:.2f}) print(f矩估计结果: μ{mu_moment:.2f}, σ{sigma_moment:.2f}) print(fMLE结果: μ{mu_mle:.2f}, σ{sigma_mle:.2f})在实际项目中我发现当样本量较大1000时两种方法结果相近但对于小样本MLE通常更可靠。一个实用的建议是先用矩估计快速获取初始值再用MLE进行精细调整。6. 性能优化与工程实践处理大规模数据时直接使用Python循环计算PDF/CDF会非常低效。我们可以利用NumPy的向量化运算进行优化def vectorized_folded_pdf(x, mu, sigma): 向量化实现的PDF计算 x np.asarray(x) mask x 0 result np.zeros_like(x) x_valid x[mask] term1 norm.pdf(x_valid, mu, sigma) term2 norm.pdf(-x_valid, mu, sigma) result[mask] term1 term2 return result def vectorized_folded_cdf(x, mu, sigma): 向量化实现的CDF计算 x np.asarray(x) mask x 0 result np.zeros_like(x) x_valid x[mask] term1 norm.cdf(x_valid, mu, sigma) term2 norm.cdf(-x_valid, mu, sigma) result[mask] term1 - term2 return result测试性能提升# 生成测试数据 large_x np.linspace(0, 10, 1000000) # 原始方法计时 %timeit [folded_normal_pdf(x, 1, 1) for x in large_x] # 向量化方法计时 %timeit vectorized_folded_pdf(large_x, 1, 1)在我的测试中向量化实现比循环快约100倍。对于需要反复计算PDF/CDF的应用如MCMC采样这种优化至关重要。另一个工程实践中的常见需求是将折叠正态分布集成到现有的统计建模流程中。我们可以创建一个符合scipy.stats接口的分布类from scipy.stats import rv_continuous class folded_normal_gen(rv_continuous): 自定义折叠正态分布类 def _pdf(self, x, mu, sigma): return vectorized_folded_pdf(x, mu, sigma) def _cdf(self, x, mu, sigma): return vectorized_folded_cdf(x, mu, sigma) def _argcheck(self, mu, sigma): return (sigma 0) folded_normal folded_normal_gen(a0, namefolded_normal)这样就能像使用标准scipy分布一样使用我们的折叠正态分布# 创建分布实例 fn folded_normal(mu1, sigma0.5) # 计算统计量 print(f均值: {fn.mean()}) print(f中位数: {fn.median()}) print(f标准差: {fn.std()}) # 生成随机数 samples fn.rvs(size1000) # 拟合数据 fit_mu, fit_sigma, _ folded_normal.fit(test_data)在最近的一个传感器数据分析项目中这种封装方式让我能够无缝地将折叠正态分布集成到现有的PyMC3贝叶斯建模流程中大大简化了开发工作。