# Python实战：数学期望与条件期望的代码实现与工程应用 在数据科学和机器学习领域，概率论中的数学期望和条件期望是构建算法模型的核心数学工具。本文将从工程实践角度出发，通过Python代码示例展示如何在实际项目中应用这些概念，帮助开发者快速掌握概率计算的编程实现。 ## 1. 数学期望的基础实现 数学期望是描述随机变量平均取值的重要指标。我们先从最基础的离散型和连续型随机变量入手，看看如何用Python进行计算。 ### 1.1 离散型随机变量的期望 对于离散型随机变量，数学期望的计算公式为： ```python def discrete_expectation(values, probabilities): """ 计算离散型随机变量的数学期望 参数: values: 随机变量取值列表 probabilities: 对应概率列表 返回: 数学期望值 """ if len(values) != len(probabilities): raise ValueError("取值和概率列表长度必须相同") if not np.isclose(sum(probabilities), 1.0): raise ValueError("概率总和必须为1") return sum(v*p for v,p in zip(values, probabilities)) # 示例：掷骰子的期望 dice_values = [1, 2, 3, 4, 5, 6] dice_probs = [1/6]*6 print(f"骰子的期望值: {discrete_expectation(dice_values, dice_probs):.2f}") ``` > 注意：实际工程中建议使用NumPy或SciPy的统计函数，这里展示原理性实现 ### 1.2 连续型随机变量的期望 连续型随机变量的期望需要通过积分计算： ```python from scipy import integrate def continuous_expectation(pdf, lower, upper): """ 计算连续型随机变量的数学期望 参数: pdf: 概率密度函数 lower: 积分下限 upper: 积分上限 返回: 数学期望值 """ # 定义被积函数：x * pdf(x) integrand = lambda x: x * pdf(x) result, _ = integrate.quad(integrand, lower, upper) return result # 示例：标准正态分布的期望 from scipy.stats import norm print(f"标准正态分布的期望: {continuous_expectation(norm.pdf, -np.inf, np.inf):.2f}") ``` ### 1.3 随机变量函数的期望 实际应用中常需要计算随机变量函数的期望： ```python def function_expectation(func, values=None, probabilities=None, pdf=None, lower=None, upper=None): """ 计算随机变量函数的期望 参数: func: 变换函数 values/probabilities: 离散情况 pdf/lower/upper: 连续情况 返回: 函数期望值 """ if values is not None: # 离散情况 transformed = [func(v) for v in values] return discrete_expectation(transformed, probabilities) elif pdf is not None: # 连续情况 integrand = lambda x: func(x) * pdf(x) result, _ = integrate.quad(integrand, lower, upper) return result else: raise ValueError("必须提供离散或连续分布参数") # 示例：计算X^2的期望，X为标准正态分布 print(f"X^2的期望(X~N(0,1)): {function_expectation(lambda x: x**2, pdf=norm.pdf, lower=-np.inf, upper=np.inf):.2f}") ``` ## 2. 条件期望的工程实现 条件期望在时间序列分析、强化学习等领域有广泛应用。下面我们看几种典型场景的实现。 ### 2.1 离散条件下的期望 当条件事件为离散情况时： ```python def conditional_expectation_discrete(values, cond_probs, given_event=None): """ 计算离散条件下的期望 参数: values: 随机变量取值 cond_probs: 条件概率字典或矩阵 given_event: 给定条件事件(可选) 返回: 条件期望值或条件期望函数 """ if given_event is not None: # 计算特定条件下的期望 probs = cond_probs[given_event] return sum(v*p for v,p in zip(values, probs)) else: # 返回条件期望函数 return lambda event: conditional_expectation_discrete(values, cond_probs, event) # 示例：天气对销售额的影响 weather_conditions = ['晴天', '雨天', '阴天'] sales_values = [100, 50, 80] cond_probs = { '晴天': [0.6, 0.3, 0.1], '雨天': [0.1, 0.7, 0.2], '阴天': [0.3, 0.4, 0.3] } print(f"雨天期望销售额: {conditional_expectation_discrete(sales_values, cond_probs, '雨天')}") ``` ### 2.2 连续条件下的期望 对于连续随机变量的条件期望，我们通常需要知道联合分布： ```python def conditional_expectation_continuous(joint_pdf, y_value, x_bounds=(-np.inf, np.inf), y_bounds=(-np.inf, np.inf)): """ 计算连续随机变量的条件期望E[X|Y=y] 参数: joint_pdf: 联合概率密度函数 y_value: Y的取值 x_bounds: X的积分范围 y_bounds: Y的积分范围(用于归一化) 返回: 条件期望值 """ # 计算条件概率密度f_X|Y(x|y) = f_XY(x,y)/f_Y(y) def conditional_pdf(x): marginal_y, _ = integrate.quad(lambda y: joint_pdf(x, y), *y_bounds) return joint_pdf(x, y_value) / marginal_y # 计算E[X|Y=y] = ∫x*f_X|Y(x|y)dx integrand = lambda x: x * conditional_pdf(x) result, _ = integrate.quad(integrand, *x_bounds) return result # 示例：二元正态分布的条件期望 from scipy.stats import multivariate_normal # 定义二元正态分布 mean = [0, 0] cov = [[1, 0.5], [0.5, 1]] rv = multivariate_normal(mean, cov) # 计算E[X|Y=0.5] print(f"E[X|Y=0.5]: {conditional_expectation_continuous(rv.pdf, 0.5):.2f}") ``` ## 3. 实际应用案例 ### 3.1 投资组合预期收益计算 假设有三种投资产品，其收益与市场状况相关： ```python # 定义市场状态和概率 market_states = ['牛市', '震荡', '熊市'] state_probs = [0.3, 0.5, 0.2] # 定义各产品在不同市场下的收益率 returns = { '股票': [0.15, 0.05, -0.1], '债券': [0.05, 0.03, 0.01], '黄金': [0.03, 0.08, 0.12] } # 计算各产品期望收益 for asset, ret in returns.items(): exp_return = discrete_expectation(ret, state_probs) print(f"{asset}期望收益率: {exp_return:.2%}") # 计算条件期望：在熊市下的期望收益 print("\n熊市条件下:") for asset, ret in returns.items(): cond_exp = discrete_expectation(ret, [0, 0, 1]) # 熊市概率为1 print(f"{asset}期望收益率: {cond_exp:.2%}") ``` ### 3.2 贝叶斯网络中的条件期望 贝叶斯网络常用于表示变量间的条件依赖关系： ```python from pgmpy.models import BayesianNetwork from pgmpy.factors.discrete import TabularCPD from pgmpy.inference import VariableElimination # 创建简单贝叶斯网络：吸烟->肺癌 model = BayesianNetwork([('Smoking', 'LungCancer')]) # 定义条件概率分布 cpd_smoking = TabularCPD('Smoking', 2, [[0.8], [0.2]]) cpd_cancer = TabularCPD('LungCancer', 2, [[0.99, 0.7], [0.01, 0.3]], evidence=['Smoking'], evidence_card=[2]) model.add_cpds(cpd_smoking, cpd_cancer) # 创建推理引擎 infer = VariableElimination(model) # 计算条件期望：P(肺癌|吸烟) q = infer.query(['LungCancer'], evidence={'Smoking': 1}) print(f"吸烟者患癌概率: {q.values[1]:.2%}") # 计算边缘期望：P(肺癌) q = infer.query(['LungCancer']) print(f"总体患癌概率: {q.values[1]:.2%}") ``` ## 4. 性能优化与工程实践 在实际工程中，我们还需要考虑计算效率和数值稳定性问题。 ### 4.1 蒙特卡洛近似方法 当解析解难以计算时，可以使用蒙特卡洛方法： ```python def monte_carlo_expectation(samples, func=None): """ 蒙特卡洛方法计算期望 参数: samples: 样本数组 func: 变换函数(可选) 返回: 期望估计值 """ if func is not None: samples = np.array([func(x) for x in samples]) return np.mean(samples) # 生成服从特定分布的样本 np.random.seed(42) samples = np.random.normal(0, 1, 10000) # 计算E[X^2] print(f"蒙特卡洛估计E[X^2]: {monte_carlo_expectation(samples, lambda x: x**2):.2f}") # 条件期望的蒙特卡洛估计 def conditional_mc(samples_x, samples_y, func, y_value, epsilon=0.1): """ 蒙特卡洛条件期望估计 参数: samples_x: X的样本 samples_y: Y的样本 func: 目标函数 y_value: 条件值 epsilon: 邻域半径 返回: 条件期望估计 """ mask = (samples_y >= y_value - epsilon) & (samples_y <= y_value + epsilon) selected = samples_x[mask] if len(selected) == 0: return np.nan return monte_carlo_expectation(selected, func) # 示例：估计E[X|Y≈0] samples_y = np.random.normal(0, 1, 10000) print(f"E[X|Y≈0]: {conditional_mc(samples, samples_y, lambda x: x, 0):.2f}") ``` ### 4.2 使用JIT加速计算 对于复杂计算，可以使用Numba进行即时编译加速： ```python from numba import jit @jit(nopython=True) def fast_discrete_expectation(values, probabilities): expectation = 0.0 for v, p in zip(values, probabilities): expectation += v * p return expectation # 大数据量测试 large_values = np.random.rand(1000000) large_probs = np.random.rand(1000000) large_probs /= large_probs.sum() # 归一化 %timeit fast_discrete_expectation(large_values, large_probs) # 通常比纯Python快10-100倍 ``` ## 5. 高级应用：强化学习中的价值函数 条件期望在强化学习的价值函数计算中扮演核心角色。我们来看一个简单的Q-learning示例： ```python import numpy as np class QLearningAgent: def __init__(self, n_states, n_actions, alpha=0.1, gamma=0.99): self.q_table = np.zeros((n_states, n_actions)) self.alpha = alpha # 学习率 self.gamma = gamma # 折扣因子 def update(self, state, action, reward, next_state, done): # 当前Q值 current_q = self.q_table[state, action] # 计算目标Q值（条件期望的估计） if done: target = reward else: max_next_q = np.max(self.q_table[next_state]) target = reward + self.gamma * max_next_q # Q值更新 self.q_table[state, action] += self.alpha * (target - current_q) def get_action(self, state, epsilon=0.1): if np.random.random() < epsilon: return np.random.randint(0, self.q_table.shape[1]) return np.argmax(self.q_table[state]) # 简单迷宫环境示例 n_states = 16 # 4x4网格 n_actions = 4 # 上下左右 agent = QLearningAgent(n_states, n_actions) # 模拟训练过程 for episode in range(1000): state = 0 # 起始状态 done = False while not done: action = agent.get_action(state) next_state, reward, done = env_step(state, action) # 假设有环境函数 agent.update(state, action, reward, next_state, done) state = next_state ``` 在这个实现中，`max_next_q`实际上是在计算当前策略下的条件期望，是强化学习中价值迭代的核心。 ## 6. 数值稳定性与异常处理 在实际工程实现中，我们需要特别注意数值稳定性问题： ```python def safe_expectation(values, probabilities=None, log_probs=None, eps=1e-10): """ 数值稳定的期望计算 参数: values: 随机变量取值 probabilities: 普通概率(可选) log_probs: 对数概率(可选，更稳定) eps: 小常数避免除零 返回: 期望值 """ if log_probs is not None: # 使用log-sum-exp技巧提高数值稳定性 max_log = np.max(log_probs) shifted_probs = np.exp(log_probs - max_log) norm = np.sum(shifted_probs) normalized_probs = shifted_probs / (norm + eps) return np.sum(values * normalized_probs) elif probabilities is not None: # 普通概率计算 probs = np.array(probabilities) probs = probs / (np.sum(probs) + eps) # 确保归一化 return np.sum(values * probs) else: raise ValueError("必须提供概率或对数概率") # 示例：处理极小概率情况 tiny_probs = [1e-30, 1e-30, 1.0-2e-30] log_probs = np.log(tiny_probs) values = [1, 2, 3] print(f"普通计算: {safe_expectation(values, probabilities=tiny_probs):.2f}") print(f"对数域计算: {safe_expectation(values, log_probs=log_probs):.2f}") ``` ## 7. 可视化分析 理解期望和条件期望的行为，可视化是非常有效的手段： ```python import matplotlib.pyplot as plt import seaborn as sns # 条件期望可视化示例 def plot_conditional_expectation(joint_pdf, y_values, x_bounds=(-3, 3)): """ 绘制条件期望E[X|Y=y]随y变化的曲线 """ expectations = [] for y in y_values: try: exp = conditional_expectation_continuous(joint_pdf, y, x_bounds) expectations.append(exp) except: expectations.append(np.nan) plt.figure(figsize=(10, 6)) plt.plot(y_values, expectations, 'b-', linewidth=2) plt.xlabel('Y value') plt.ylabel('E[X|Y=y]') plt.title('Conditional Expectation Function') plt.grid(True) plt.show() # 定义联合分布(二元正态) def bivariate_normal_pdf(x, y): return multivariate_normal([0, 0], [[1, 0.8], [0.8, 1]]).pdf([x, y]) # 绘制条件期望曲线 y_range = np.linspace(-2, 2, 50) plot_conditional_expectation(bivariate_normal_pdf, y_range) ``` ## 8. 测试与验证 确保我们的实现正确性至关重要： ```python import unittest class TestExpectations(unittest.TestCase): def setUp(self): self.values = [1, 2, 3, 4] self.probs = [0.1, 0.2, 0.3, 0.4] self.pdf = norm(0, 1).pdf self.samples = norm(0, 1).rvs(10000, random_state=42) def test_discrete_expectation(self): expected = sum(v*p for v,p in zip(self.values, self.probs)) result = discrete_expectation(self.values, self.probs) self.assertAlmostEqual(expected, result, places=6) def test_continuous_expectation(self): expected = 0 # 标准正态期望 result = continuous_expectation(self.pdf, -np.inf, np.inf) self.assertAlmostEqual(expected, result, places=2) def test_monte_carlo(self): result = monte_carlo_expectation(self.samples) self.assertTrue(-0.1 < result < 0.1) # 应该在0附近 if __name__ == '__main__': unittest.main(argv=[''], exit=False) ``` ## 9. 实际工程建议 在真实项目中应用这些概念时，有几个实用建议： 1. **优先使用成熟库函数**：对于常见分布，优先使用SciPy或NumPy的内置函数 ```python from scipy.stats import norm mean, var = norm.stats(moments='mv') ``` 2. **处理边界情况**：特别是连续分布，注意积分限和数值稳定性 3. **缓存中间结果**：对于昂贵的条件期望计算，考虑缓存结果 4. **并行化计算**：对于蒙特卡洛等可并行算法，使用多进程加速 ```python from multiprocessing import Pool def parallel_mc(samples, func, n_workers=4): with Pool(n_workers) as p: results = p.map(func, samples) return np.mean(results) ``` 5. **考虑对数空间计算**：概率相乘容易下溢，转为对数空间相加更稳定 ## 10. 扩展应用：金融衍生品定价 在金融工程中，条件期望用于衍生品定价。以欧式看涨期权为例： ```python def european_call_price(S0, K, T, r, sigma, n_simulations=100000): """ 蒙特卡洛方法计算欧式看涨期权价格 参数: S0: 标的资产现价 K: 行权价 T: 到期时间(年) r: 无风险利率 sigma: 波动率 n_simulations: 模拟次数 返回: 期权价格 """ # 生成标的资产价格路径 z = np.random.normal(0, 1, n_simulations) ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*z) # 计算payoff的条件期望并贴现 payoffs = np.maximum(ST - K, 0) return np.exp(-r*T) * np.mean(payoffs) # 示例计算 price = european_call_price(S0=100, K=105, T=1, r=0.05, sigma=0.2) print(f"欧式看涨期权价格: {price:.2f}") ``` 这个例子中，期权价格实际上是未来收益的条件期望的现值，展示了条件期望在金融工程中的核心作用。