pengxiao
219232de74
以讲义内容为骨架迁移到标准目录格式： - officefile/ 主内容（12章 + 附录 + CC4SI补充） - dofile/ 代码示例（11个Python脚本） - data/ 图片资源 - output/ 生成输出（忽略） - Archive/ 归档旧目录（忽略） - .claude/skills/ 保留markdown-to-docx工具链 - .pandoc/ 保留CSL和本地化配置 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
939 lines
29 KiB
Python
939 lines
29 KiB
Python
"""
|
||
概率与不确定性示例 (Probability and Uncertainty Example)
|
||
========================================================
|
||
|
||
本示例展示如何在空间智能系统中处理概率和不确定性。
|
||
在空间决策中，不确定性是普遍存在的，理解和管理不确定性
|
||
对于做出可靠的决策至关重要。
|
||
|
||
核心概念:
|
||
1. 概率分布 - 描述随机变量的可能取值及其概率
|
||
2. 贝叶斯推理 - 基于新证据更新信念
|
||
3. 蒙特卡洛模拟 - 通过随机采样评估不确定性
|
||
4. 置信区间 - 估计结果的范围
|
||
5. 敏感性分析 - 评估输入变化对输出的影响
|
||
|
||
应用场景:
|
||
- 空间插值的不确定性量化
|
||
- 多准则决策的敏感性分析
|
||
- 风险评估与概率预测
|
||
- 传感器数据的可靠性分析
|
||
|
||
作者: CC4SI 项目组
|
||
"""
|
||
|
||
import math
|
||
import random
|
||
from typing import List, Dict, Tuple, Optional, Callable, Any
|
||
from dataclasses import dataclass, field
|
||
from abc import ABC, abstractmethod
|
||
from enum import Enum
|
||
import statistics
|
||
|
||
|
||
# ============================================================================
|
||
# 概率分布基础类
|
||
# ============================================================================
|
||
|
||
class DistributionType(Enum):
|
||
"""分布类型枚举"""
|
||
NORMAL = "normal" # 正态分布
|
||
UNIFORM = "uniform" # 均匀分布
|
||
TRIANGULAR = "triangular" # 三角分布
|
||
EXPONENTIAL = "exponential" # 指数分布
|
||
BETA = "beta" # Beta分布
|
||
|
||
|
||
class ProbabilityDistribution(ABC):
|
||
"""概率分布抽象基类"""
|
||
|
||
def __init__(self, name: str = ""):
|
||
self.name = name
|
||
|
||
@abstractmethod
|
||
def sample(self) -> float:
|
||
"""从分布中采样一个值"""
|
||
pass
|
||
|
||
@abstractmethod
|
||
def mean(self) -> float:
|
||
"""计算期望值"""
|
||
pass
|
||
|
||
@abstractmethod
|
||
def std(self) -> float:
|
||
"""计算标准差"""
|
||
pass
|
||
|
||
@abstractmethod
|
||
def pdf(self, x: float) -> float:
|
||
"""概率密度函数"""
|
||
pass
|
||
|
||
def cdf(self, x: float) -> float:
|
||
"""累积分布函数 (近似计算)"""
|
||
# 使用蒙特卡洛积分近似
|
||
n_samples = 10000
|
||
count = sum(1 for _ in range(n_samples) if self.sample() <= x)
|
||
return count / n_samples
|
||
|
||
def confidence_interval(self, confidence: float = 0.95) -> Tuple[float, float]:
|
||
"""计算置信区间"""
|
||
n_samples = 10000
|
||
samples = [self.sample() for _ in range(n_samples)]
|
||
alpha = 1 - confidence
|
||
lower = quantile(samples, alpha / 2)
|
||
upper = quantile(samples, 1 - alpha / 2)
|
||
return lower, upper
|
||
|
||
|
||
def quantile(data: List[float], q: float) -> float:
|
||
"""计算分位数"""
|
||
sorted_data = sorted(data)
|
||
index = int(q * len(sorted_data))
|
||
return sorted_data[min(index, len(sorted_data) - 1)]
|
||
|
||
|
||
# ============================================================================
|
||
# 具体概率分布实现
|
||
# ============================================================================
|
||
|
||
class NormalDistribution(ProbabilityDistribution):
|
||
"""
|
||
正态分布 (高斯分布)
|
||
|
||
最常用的连续概率分布，由均值和标准差参数化。
|
||
"""
|
||
|
||
def __init__(self, mu: float = 0.0, sigma: float = 1.0, name: str = ""):
|
||
super().__init__(name)
|
||
self.mu = mu # 均值
|
||
self.sigma = sigma # 标准差
|
||
if sigma <= 0:
|
||
raise ValueError("标准差必须为正数")
|
||
|
||
def sample(self) -> float:
|
||
"""使用 Box-Muller 变换生成正态分布随机数"""
|
||
u1 = random.random()
|
||
u2 = random.random()
|
||
while u1 == 0: # 避免log(0)
|
||
u1 = random.random()
|
||
z0 = math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)
|
||
return self.mu + self.sigma * z0
|
||
|
||
def mean(self) -> float:
|
||
return self.mu
|
||
|
||
def std(self) -> float:
|
||
return self.sigma
|
||
|
||
def pdf(self, x: float) -> float:
|
||
"""正态分布概率密度函数"""
|
||
coeff = 1 / (self.sigma * math.sqrt(2 * math.pi))
|
||
exponent = -0.5 * ((x - self.mu) / self.sigma) ** 2
|
||
return coeff * math.exp(exponent)
|
||
|
||
def __repr__(self) -> str:
|
||
return f"Normal(μ={self.mu}, σ={self.sigma})"
|
||
|
||
|
||
class UniformDistribution(ProbabilityDistribution):
|
||
"""
|
||
均匀分布
|
||
|
||
在指定范围内等概率取值。
|
||
"""
|
||
|
||
def __init__(self, a: float = 0.0, b: float = 1.0, name: str = ""):
|
||
super().__init__(name)
|
||
self.a = a # 下界
|
||
self.b = b # 上界
|
||
if a >= b:
|
||
raise ValueError("下界必须小于上界")
|
||
|
||
def sample(self) -> float:
|
||
return self.a + (self.b - self.a) * random.random()
|
||
|
||
def mean(self) -> float:
|
||
return (self.a + self.b) / 2
|
||
|
||
def std(self) -> float:
|
||
return (self.b - self.a) / math.sqrt(12)
|
||
|
||
def pdf(self, x: float) -> float:
|
||
if self.a <= x <= self.b:
|
||
return 1 / (self.b - self.a)
|
||
return 0.0
|
||
|
||
def __repr__(self) -> str:
|
||
return f"Uniform({self.a}, {self.b})"
|
||
|
||
|
||
class TriangularDistribution(ProbabilityDistribution):
|
||
"""
|
||
三角分布
|
||
|
||
由最小值、最大值和众数定义的分布，常用于
|
||
当只知道边界和最可能值时建模不确定性。
|
||
"""
|
||
|
||
def __init__(self, a: float, b: float, c: float, name: str = ""):
|
||
super().__init__(name)
|
||
self.a = a # 最小值
|
||
self.b = b # 最大值
|
||
self.c = c # 众数 (最可能值)
|
||
if not (a <= c <= b):
|
||
raise ValueError("必须满足 a <= c <= b")
|
||
|
||
def sample(self) -> float:
|
||
u = random.random()
|
||
fc = (self.c - self.a) / (self.b - self.a)
|
||
if u < fc:
|
||
return self.a + math.sqrt(u * (self.b - self.a) * (self.c - self.a))
|
||
else:
|
||
return self.b - math.sqrt((1 - u) * (self.b - self.a) * (self.b - self.c))
|
||
|
||
def mean(self) -> float:
|
||
return (self.a + self.b + self.c) / 3
|
||
|
||
def std(self) -> float:
|
||
numerator = (self.a**2 + self.b**2 + self.c**2
|
||
- self.a * self.b - self.a * self.c - self.b * self.c)
|
||
return math.sqrt(numerator / 18)
|
||
|
||
def pdf(self, x: float) -> float:
|
||
if x < self.a or x > self.b:
|
||
return 0.0
|
||
if x < self.c:
|
||
return 2 * (x - self.a) / ((self.b - self.a) * (self.c - self.a))
|
||
else:
|
||
return 2 * (self.b - x) / ((self.b - self.a) * (self.b - self.c))
|
||
|
||
def __repr__(self) -> str:
|
||
return f"Triangular({self.a}, {self.c}, {self.b})"
|
||
|
||
|
||
# ============================================================================
|
||
# 贝叶斯推理
|
||
# ============================================================================
|
||
|
||
@dataclass
|
||
class BayesianBelief:
|
||
"""
|
||
贝叶斯信念状态
|
||
|
||
表示对某个假设的信念，包含先验、似然和后验。
|
||
"""
|
||
hypothesis: str
|
||
prior: float # 先验概率 P(H)
|
||
likelihood: float # 似然 P(E|H)
|
||
evidence: Optional[float] = None # 证据概率 P(E)
|
||
posterior: Optional[float] = None # 后验概率 P(H|E)
|
||
|
||
def update(self, evidence_prob: float = None) -> float:
|
||
"""
|
||
更新后验概率
|
||
|
||
Args:
|
||
evidence_prob: P(E)，如果None则使用归一化
|
||
|
||
Returns:
|
||
后验概率
|
||
"""
|
||
# P(H|E) = P(E|H) * P(H) / P(E)
|
||
numerator = self.likelihood * self.prior
|
||
|
||
if evidence_prob is not None:
|
||
self.evidence = evidence_prob
|
||
self.posterior = numerator / evidence_prob
|
||
else:
|
||
# 假设有多个假设，需要归一化
|
||
self.posterior = numerator # 简化版本
|
||
|
||
return self.posterior
|
||
|
||
|
||
class BayesianUpdater:
|
||
"""
|
||
贝叶斯更新器
|
||
|
||
管理多个假设的贝叶斯更新。
|
||
"""
|
||
|
||
def __init__(self, hypotheses: List[str]):
|
||
"""
|
||
初始化贝叶斯更新器
|
||
|
||
Args:
|
||
hypotheses: 假设列表
|
||
"""
|
||
# 初始化先验概率 (均匀分布)
|
||
prior = 1.0 / len(hypotheses)
|
||
self.beliefs: Dict[str, BayesianBelief] = {
|
||
h: BayesianBelief(hypothesis=h, prior=prior, likelihood=1.0)
|
||
for h in hypotheses
|
||
}
|
||
|
||
def set_prior(self, hypothesis: str, prior: float) -> None:
|
||
"""设置先验概率"""
|
||
if hypothesis in self.beliefs:
|
||
self.beliefs[hypothesis].prior = prior
|
||
|
||
def update_with_evidence(self, likelihoods: Dict[str, float]) -> None:
|
||
"""
|
||
用证据更新所有假设
|
||
|
||
Args:
|
||
likelihoods: 每个假设的似然 P(E|H)
|
||
"""
|
||
# 更新似然
|
||
for h, likelihood in likelihoods.items():
|
||
if h in self.beliefs:
|
||
self.beliefs[h].likelihood = likelihood
|
||
|
||
# 计算证据概率 (归一化常数)
|
||
evidence = sum(
|
||
b.likelihood * b.prior
|
||
for b in self.beliefs.values()
|
||
)
|
||
|
||
# 更新后验
|
||
for belief in self.beliefs.values():
|
||
belief.update(evidence)
|
||
|
||
def get_posteriors(self) -> Dict[str, float]:
|
||
"""获取所有后验概率"""
|
||
return {
|
||
h: b.posterior or b.prior
|
||
for h, b in self.beliefs.items()
|
||
}
|
||
|
||
def get_most_likely(self) -> Tuple[str, float]:
|
||
"""获取最可能的假设"""
|
||
posteriors = self.get_posteriors()
|
||
return max(posteriors.items(), key=lambda x: x[1])
|
||
|
||
def print_beliefs(self) -> None:
|
||
"""打印信念状态"""
|
||
print("\n贝叶斯信念状态:")
|
||
print("-" * 60)
|
||
print(f"{'假设':<20} {'先验':<12} {'似然':<12} {'后验':<12}")
|
||
print("-" * 60)
|
||
for belief in self.beliefs.values():
|
||
posterior = belief.posterior if belief.posterior is not None else belief.prior
|
||
print(f"{belief.hypothesis:<20} {belief.prior:<12.4f} "
|
||
f"{belief.likelihood:<12.4f} {posterior:<12.4f}")
|
||
print("-" * 60)
|
||
|
||
|
||
# ============================================================================
|
||
# 蒙特卡洛模拟
|
||
# ============================================================================
|
||
|
||
@dataclass
|
||
class SimulationResult:
|
||
"""模拟结果"""
|
||
samples: List[float] = field(default_factory=list)
|
||
mean: float = 0.0
|
||
std: float = 0.0
|
||
min: float = 0.0
|
||
max: float = 0.0
|
||
median: float = 0.0
|
||
confidence_interval: Tuple[float, float] = (0.0, 0.0)
|
||
percentiles: Dict[float, float] = field(default_factory=dict)
|
||
|
||
def calculate_statistics(self, confidence: float = 0.95) -> None:
|
||
"""计算统计量"""
|
||
if not self.samples:
|
||
return
|
||
|
||
self.mean = statistics.mean(self.samples)
|
||
self.std = statistics.stdev(self.samples) if len(self.samples) > 1 else 0
|
||
self.min = min(self.samples)
|
||
self.max = max(self.samples)
|
||
self.median = statistics.median(self.samples)
|
||
|
||
# 置信区间
|
||
alpha = 1 - confidence
|
||
sorted_samples = sorted(self.samples)
|
||
n = len(sorted_samples)
|
||
self.confidence_interval = (
|
||
sorted_samples[int(alpha / 2 * n)],
|
||
sorted_samples[int((1 - alpha / 2) * n)]
|
||
)
|
||
|
||
# 常用百分位数
|
||
for p in [0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]:
|
||
self.percentiles[p] = sorted_samples[int(p * n)]
|
||
|
||
def print_summary(self) -> None:
|
||
"""打印结果摘要"""
|
||
print(f"\n蒙特卡洛模拟结果 (n={len(self.samples)}):")
|
||
print("-" * 50)
|
||
print(f"均值: {self.mean:.4f}")
|
||
print(f"中位数: {self.median:.4f}")
|
||
print(f"标准差: {self.std:.4f}")
|
||
print(f"范围: [{self.min:.4f}, {self.max:.4f}]")
|
||
print(f"95% 置信区间: [{self.confidence_interval[0]:.4f}, "
|
||
f"{self.confidence_interval[1]:.4f}]")
|
||
print(f"\n百分位数:")
|
||
for p, value in sorted(self.percentiles.items()):
|
||
print(f" {p*100:>5.0f}%: {value:.4f}")
|
||
print("-" * 50)
|
||
|
||
|
||
class MonteCarloSimulator:
|
||
"""
|
||
蒙特卡洛模拟器
|
||
|
||
通过随机采样评估不确定性。
|
||
"""
|
||
|
||
def __init__(self, seed: int = None):
|
||
"""初始化模拟器"""
|
||
if seed is not None:
|
||
random.seed(seed)
|
||
|
||
def simulate(self, model: Callable[[], float],
|
||
n_runs: int = 10000) -> SimulationResult:
|
||
"""
|
||
运行模拟
|
||
|
||
Args:
|
||
model: 返回模拟值的函数
|
||
n_runs: 运行次数
|
||
|
||
Returns:
|
||
模拟结果
|
||
"""
|
||
samples = [model() for _ in range(n_runs)]
|
||
result = SimulationResult(samples=samples)
|
||
result.calculate_statistics()
|
||
return result
|
||
|
||
def simulate_with_inputs(self,
|
||
model: Callable[[Dict[str, float]], float],
|
||
input_distributions: Dict[str, ProbabilityDistribution],
|
||
n_runs: int = 10000) -> SimulationResult:
|
||
"""
|
||
使用输入分布运行模拟
|
||
|
||
Args:
|
||
model: 接受输入字典的函数
|
||
input_distributions: 输入变量到其分布的映射
|
||
n_runs: 运行次数
|
||
|
||
Returns:
|
||
模拟结果
|
||
"""
|
||
samples = []
|
||
for _ in range(n_runs):
|
||
inputs = {
|
||
name: dist.sample()
|
||
for name, dist in input_distributions.items()
|
||
}
|
||
samples.append(model(inputs))
|
||
|
||
result = SimulationResult(samples=samples)
|
||
result.calculate_statistics()
|
||
return result
|
||
|
||
|
||
# ============================================================================
|
||
# 敏感性分析
|
||
# ============================================================================
|
||
|
||
@dataclass
|
||
class SensitivityResult:
|
||
"""敏感性分析结果"""
|
||
sensitivity_coefficients: Dict[str, float] = field(default_factory=dict)
|
||
rankings: List[Tuple[str, float]] = field(default_factory=list)
|
||
tornado_data: Dict[str, Tuple[float, float]] = field(default_factory=dict)
|
||
|
||
def print_summary(self) -> None:
|
||
"""打印敏感性分析摘要"""
|
||
print("\n敏感性分析结果:")
|
||
print("-" * 50)
|
||
print("排名 | 变量 | 敏感性系数")
|
||
print("-" * 50)
|
||
for i, (var, coef) in enumerate(self.rankings, 1):
|
||
print(f"{i:4d} | {var:<11} | {coef:10.4f}")
|
||
print("-" * 50)
|
||
|
||
|
||
class SensitivityAnalyzer:
|
||
"""
|
||
敏感性分析器
|
||
|
||
评估输入变化对输出的影响。
|
||
"""
|
||
|
||
def __init__(self):
|
||
self.model: Optional[Callable] = None
|
||
self.base_inputs: Optional[Dict[str, float]] = None
|
||
|
||
def simple_sensitivity(self,
|
||
model: Callable[[Dict[str, float]], float],
|
||
base_inputs: Dict[str, float],
|
||
variations: Dict[str, float] = None) -> SensitivityResult:
|
||
"""
|
||
简单敏感性分析 (单因素)
|
||
|
||
Args:
|
||
model: 待分析模型
|
||
base_inputs: 基准输入值
|
||
variations: 各变量的变化幅度 (默认 ±10%)
|
||
|
||
Returns:
|
||
敏感性结果
|
||
"""
|
||
if variations is None:
|
||
variations = {k: 0.1 for k in base_inputs.keys()}
|
||
|
||
# 计算基准输出
|
||
base_output = model(base_inputs)
|
||
|
||
# 计算敏感性系数
|
||
coefficients = {}
|
||
tornado_data = {}
|
||
|
||
for var, variation in variations.items():
|
||
original_value = base_inputs[var]
|
||
|
||
# 正向变化
|
||
base_inputs[var] = original_value * (1 + variation)
|
||
output_plus = model(base_inputs)
|
||
|
||
# 负向变化
|
||
base_inputs[var] = original_value * (1 - variation)
|
||
output_minus = model(base_inputs)
|
||
|
||
# 恢复原值
|
||
base_inputs[var] = original_value
|
||
|
||
# 计算敏感性系数 (归一化)
|
||
delta_output = output_plus - output_minus
|
||
delta_input = 2 * variation * original_value
|
||
coefficient = delta_output / delta_input if delta_input != 0 else 0
|
||
|
||
coefficients[var] = coefficient
|
||
tornado_data[var] = (output_minus, output_plus)
|
||
|
||
# 排名
|
||
rankings = sorted(coefficients.items(), key=lambda x: abs(x[1]), reverse=True)
|
||
|
||
return SensitivityResult(
|
||
sensitivity_coefficients=coefficients,
|
||
rankings=rankings,
|
||
tornado_data=tornado_data
|
||
)
|
||
|
||
def regression_sensitivity(self,
|
||
model: Callable[[Dict[str, float]], float],
|
||
input_distributions: Dict[str, ProbabilityDistribution],
|
||
n_samples: int = 1000) -> SensitivityResult:
|
||
"""
|
||
基于回归的敏感性分析
|
||
|
||
Args:
|
||
model: 待分析模型
|
||
input_distributions: 输入分布
|
||
n_samples: 样本数量
|
||
|
||
Returns:
|
||
敏感性结果
|
||
"""
|
||
simulator = MonteCarloSimulator()
|
||
|
||
# 生成样本
|
||
input_samples = []
|
||
output_samples = []
|
||
|
||
for _ in range(n_samples):
|
||
inputs = {
|
||
name: dist.sample()
|
||
for name, dist in input_distributions.items()
|
||
}
|
||
input_samples.append(inputs)
|
||
output_samples.append(model(inputs))
|
||
|
||
# 计算标准化回归系数 (SRC)
|
||
# SRC = beta * (std_x / std_y)
|
||
|
||
import statistics
|
||
|
||
std_y = statistics.stdev(output_samples)
|
||
coefficients = {}
|
||
|
||
for var in input_distributions.keys():
|
||
x_values = [s[var] for s in input_samples]
|
||
std_x = statistics.stdev(x_values)
|
||
|
||
# 计算相关系数
|
||
mean_x = statistics.mean(x_values)
|
||
mean_y = statistics.mean(output_samples)
|
||
|
||
numerator = sum((x - mean_x) * (y - mean_y)
|
||
for x, y in zip(x_values, output_samples))
|
||
denominator = math.sqrt(
|
||
sum((x - mean_x)**2 for x in x_values) *
|
||
sum((y - mean_y)**2 for y in output_samples)
|
||
)
|
||
|
||
correlation = numerator / denominator if denominator != 0 else 0
|
||
src = correlation * (std_x / std_y) if std_y > 0 else 0
|
||
|
||
coefficients[var] = src
|
||
|
||
rankings = sorted(coefficients.items(), key=lambda x: abs(x[1]), reverse=True)
|
||
|
||
return SensitivityResult(
|
||
sensitivity_coefficients=coefficients,
|
||
rankings=rankings
|
||
)
|
||
|
||
|
||
# ============================================================================
|
||
# 空间概率应用示例
|
||
# ============================================================================
|
||
|
||
class SpatialProbabilityModel:
|
||
"""
|
||
空间概率模型
|
||
|
||
将概率理论应用于空间问题。
|
||
"""
|
||
|
||
@staticmethod
|
||
def uncertain_distance(point1: Tuple[float, float],
|
||
point2: Tuple[float, float],
|
||
distance_error_std: float = 5.0) -> NormalDistribution:
|
||
"""
|
||
带不确定性的距离计算
|
||
|
||
Args:
|
||
point1: 第一个点 (x, y)
|
||
point2: 第二个点 (x, y)
|
||
distance_error_std: 距离测量误差的标准差
|
||
|
||
Returns:
|
||
距离的概率分布
|
||
"""
|
||
# 计算确定性距离
|
||
dx = point2[0] - point1[0]
|
||
dy = point2[1] - point1[1]
|
||
true_distance = math.sqrt(dx**2 + dy**2)
|
||
|
||
# 返回正态分布
|
||
return NormalDistribution(mu=true_distance, sigma=distance_error_std)
|
||
|
||
@staticmethod
|
||
def location_probability(measurement: Tuple[float, float],
|
||
true_location: Tuple[float, float],
|
||
measurement_error: float = 10.0) -> float:
|
||
"""
|
||
计算测量位置的似然概率
|
||
|
||
Args:
|
||
measurement: 测量位置
|
||
true_location: 真实位置
|
||
measurement_error: 测量误差标准差
|
||
|
||
Returns:
|
||
似然概率
|
||
"""
|
||
dist = SpatialProbabilityModel.uncertain_distance(
|
||
measurement, true_location, measurement_error
|
||
)
|
||
# 使用正态分布 PDF
|
||
return dist.pdf(0)
|
||
|
||
@staticmethod
|
||
def bayesian_location_update(prior_locations: List[Tuple[float, float]],
|
||
measurements: List[Tuple[float, float]],
|
||
measurement_error: float = 10.0) -> List[float]:
|
||
"""
|
||
贝叶斯位置更新
|
||
|
||
Args:
|
||
prior_locations: 候选真实位置列表
|
||
measurements: 测量位置列表
|
||
measurement_error: 测量误差
|
||
|
||
Returns:
|
||
每个候选位置的后验概率
|
||
"""
|
||
n = len(prior_locations)
|
||
posteriors = []
|
||
|
||
for candidate in prior_locations:
|
||
# 计算似然 (所有测量的乘积)
|
||
likelihood = 1.0
|
||
for measurement in measurements:
|
||
prob = SpatialProbabilityModel.location_probability(
|
||
measurement, candidate, measurement_error
|
||
)
|
||
likelihood *= prob
|
||
|
||
# 先验 (均匀)
|
||
prior = 1.0 / n
|
||
|
||
# 后验 (未归一化)
|
||
posterior = likelihood * prior
|
||
posteriors.append(posterior)
|
||
|
||
# 归一化
|
||
total = sum(posteriors)
|
||
if total > 0:
|
||
posteriors = [p / total for p in posteriors]
|
||
|
||
return posteriors
|
||
|
||
|
||
# ============================================================================
|
||
# 主程序
|
||
# ============================================================================
|
||
|
||
def main():
|
||
"""主程序 - 演示概率与不确定性的使用"""
|
||
|
||
print("="*70)
|
||
print("概率与不确定性示例演示")
|
||
print("="*70)
|
||
|
||
random.seed(42) # 可重现的结果
|
||
|
||
# ========================================================================
|
||
# 1. 概率分布示例
|
||
# ========================================================================
|
||
print("\n[部分 1] 概率分布")
|
||
print("-" * 50)
|
||
|
||
# 创建不同的分布
|
||
normal = NormalDistribution(mu=100, sigma=15, name="温度")
|
||
uniform = UniformDistribution(a=50, b=150, name="范围")
|
||
triangular = TriangularDistribution(a=60, c=100, b=140, name="估计")
|
||
|
||
distributions = [normal, uniform, triangular]
|
||
|
||
for dist in distributions:
|
||
print(f"\n{dist}:")
|
||
print(f" 均值: {dist.mean():.2f}")
|
||
print(f" 标准差: {dist.std():.2f}")
|
||
samples = [dist.sample() for _ in range(5)]
|
||
print(f" 样本: {[f'{s:.2f}' for s in samples]}")
|
||
|
||
ci = dist.confidence_interval(0.95)
|
||
print(f" 95% CI: [{ci[0]:.2f}, {ci[1]:.2f}]")
|
||
|
||
# ========================================================================
|
||
# 2. 贝叶斯推理示例
|
||
# ========================================================================
|
||
print("\n\n[部分 2] 贝叶斯推理")
|
||
print("-" * 50)
|
||
print("问题: 根据土壤测试结果判断土地适宜性")
|
||
|
||
# 假设: 土地适宜性等级
|
||
hypotheses = ["高适宜", "中适宜", "低适宜", "不适宜"]
|
||
updater = BayesianUpdater(hypotheses)
|
||
|
||
# 设置不同的先验 (基于历史数据)
|
||
updater.set_prior("高适宜", 0.2)
|
||
updater.set_prior("中适宜", 0.3)
|
||
updater.set_prior("低适宜", 0.3)
|
||
updater.set_prior("不适宜", 0.2)
|
||
|
||
print("\n初始信念:")
|
||
updater.print_beliefs()
|
||
|
||
# 证据1: 土壤pH值检测
|
||
print("\n证据1: 土壤pH值适中 (似然更新)")
|
||
updater.update_with_evidence({
|
||
"高适宜": 0.8, # pH值对高适宜的可能性高
|
||
"中适宜": 0.6,
|
||
"低适宜": 0.3,
|
||
"不适宜": 0.1
|
||
})
|
||
updater.print_beliefs()
|
||
|
||
# 证据2: 有机质含量检测
|
||
print("\n证据2: 有机质含量高 (似然更新)")
|
||
updater.update_with_evidence({
|
||
"高适宜": 0.9,
|
||
"中适宜": 0.5,
|
||
"低适宜": 0.2,
|
||
"不适宜": 0.05
|
||
})
|
||
updater.print_beliefs()
|
||
|
||
most_likely = updater.get_most_likely()
|
||
print(f"\n最可能的假设: {most_likely[0]} (概率: {most_likely[1]:.2%})")
|
||
|
||
# ========================================================================
|
||
# 3. 蒙特卡洛模拟示例
|
||
# ========================================================================
|
||
print("\n\n[部分 3] 蒙特卡洛模拟")
|
||
print("-" * 50)
|
||
print("问题: 评估房地产开发项目的预期收益")
|
||
|
||
def development_model(inputs: Dict[str, float]) -> float:
|
||
"""房地产开发收益模型"""
|
||
land_cost = inputs["land_cost"]
|
||
construction_cost = inputs["construction_cost"]
|
||
selling_price = inputs["selling_price"]
|
||
units = inputs["units"]
|
||
sales_rate = inputs["sales_rate"]
|
||
|
||
# 收益 = (售价 * 单元数 * 销售率) - (土地成本 + 建设成本 * 单元数)
|
||
revenue = selling_price * units * sales_rate
|
||
total_cost = land_cost + construction_cost * units
|
||
return revenue - total_cost
|
||
|
||
# 定义输入分布
|
||
input_dists = {
|
||
"land_cost": TriangularDistribution(800000, 1000000, 1500000), # 土地成本
|
||
"construction_cost": NormalDistribution(50000, 5000), # 单元建设成本
|
||
"selling_price": NormalDistribution(150000, 15000), # 单元售价
|
||
"units": TriangularDistribution(80, 100, 120), # 单元数量
|
||
"sales_rate": BetaDistribution(alpha=8, beta=2, a=0, b=1) # 销售率
|
||
}
|
||
|
||
simulator = MonteCarloSimulator()
|
||
result = simulator.simulate_with_inputs(development_model, input_dists, n_runs=10000)
|
||
|
||
result.print_summary()
|
||
|
||
# 风险评估
|
||
negative_prob = sum(1 for s in result.samples if s < 0) / len(result.samples)
|
||
print(f"\n风险分析:")
|
||
print(f" 亏损概率: {negative_prob:.2%}")
|
||
profit_prob = sum(1 for s in result.samples if s > 1000000) / len(result.samples)
|
||
print(f" 超过100万利润概率: {profit_prob:.2%}")
|
||
|
||
# ========================================================================
|
||
# 4. 敏感性分析示例
|
||
# ========================================================================
|
||
print("\n\n[部分 4] 敏感性分析")
|
||
print("-" * 50)
|
||
print("问题: 分析各因素对收益的影响程度")
|
||
|
||
analyzer = SensitivityAnalyzer()
|
||
|
||
# 简单敏感性分析
|
||
base_inputs = {
|
||
"land_cost": 1000000,
|
||
"construction_cost": 50000,
|
||
"selling_price": 150000,
|
||
"units": 100,
|
||
"sales_rate": 0.85
|
||
}
|
||
|
||
sensitivity_result = analyzer.simple_sensitivity(
|
||
development_model, base_inputs, variations={k: 0.1 for k in base_inputs.keys()}
|
||
)
|
||
|
||
sensitivity_result.print_summary()
|
||
|
||
# ========================================================================
|
||
# 5. 空间概率应用
|
||
# ========================================================================
|
||
print("\n\n[部分 5] 空间概率应用")
|
||
print("-" * 50)
|
||
print("问题: GPS定位的不确定性")
|
||
|
||
# 真实位置
|
||
true_location = (1000, 2000)
|
||
|
||
# 带误差的测量
|
||
measurements = [
|
||
(1005, 2003),
|
||
(998, 1998),
|
||
(1002, 2005),
|
||
(995, 2000)
|
||
]
|
||
|
||
# 候选位置
|
||
candidates = [
|
||
(1000, 2000), # 真实位置
|
||
(1015, 2015),
|
||
(990, 1990),
|
||
(1005, 1995)
|
||
]
|
||
|
||
posteriors = SpatialProbabilityModel.bayesian_location_update(
|
||
candidates, measurements, measurement_error=5.0
|
||
)
|
||
|
||
print("\n候选位置的后验概率:")
|
||
for i, (loc, prob) in enumerate(zip(candidates, posteriors)):
|
||
print(f" 位置 {i+1} {loc}: {prob:.4f}")
|
||
|
||
most_likely_idx = max(range(len(posteriors)), key=lambda i: posteriors[i])
|
||
print(f"\n最可能的位置: 位置 {most_likely_idx+1} {candidates[most_likely_idx]}")
|
||
|
||
print("\n" + "="*70)
|
||
print("演示完成!")
|
||
print("="*70)
|
||
|
||
|
||
# Beta分布实现 (用于上面代码中的引用)
|
||
class BetaDistribution(ProbabilityDistribution):
|
||
"""Beta分布 - 用于建模[0,1]区间内的概率"""
|
||
|
||
def __init__(self, alpha: float, beta: float, a: float = 0, b: float = 1, name: str = ""):
|
||
super().__init__(name)
|
||
self.alpha = alpha
|
||
self.beta = beta
|
||
self.a = a # 下界
|
||
self.b = b # 上界
|
||
|
||
def sample(self) -> float:
|
||
# 使用numpy的gamma函数近似
|
||
import math
|
||
import random
|
||
|
||
# 生成Gamma随机变量
|
||
def gamma(alpha):
|
||
if alpha < 1:
|
||
return gamma(alpha + 1) * (random.random() ** (1 / alpha))
|
||
# Marsaglia and Tsang's method
|
||
d = alpha - 1/3
|
||
c = 1 / math.sqrt(9 * d)
|
||
while True:
|
||
x = random.gauss(0, 1)
|
||
v = (1 + c * x) ** 3
|
||
if v > 0:
|
||
u = random.random()
|
||
if u < 1 - 0.0331 * (x * x) ** 2:
|
||
return d * v
|
||
if math.log(u) < 0.5 * x * x + d * (1 - v + math.log(v)):
|
||
return d * v
|
||
|
||
x = gamma(self.alpha)
|
||
y = gamma(self.beta)
|
||
beta_sample = x / (x + y)
|
||
|
||
# 转换到[a, b]区间
|
||
return self.a + (self.b - self.a) * beta_sample
|
||
|
||
def mean(self) -> float:
|
||
return self.a + (self.b - self.a) * self.alpha / (self.alpha + self.beta)
|
||
|
||
def std(self) -> float:
|
||
mean_raw = self.alpha / (self.alpha + self.beta)
|
||
var_raw = (self.alpha * self.beta) / (
|
||
(self.alpha + self.beta) ** 2 * (self.alpha + self.beta + 1)
|
||
)
|
||
return (self.b - self.a) * math.sqrt(var_raw)
|
||
|
||
def pdf(self, x: float) -> float:
|
||
# 简化版本，仅返回近似值
|
||
return 1.0 # 实际应实现Beta分布的PDF
|
||
|
||
def __repr__(self) -> str:
|
||
return f"Beta(α={self.alpha}, β={self.beta})"
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|