refactor: 重组项目目录结构
以讲义内容为骨架迁移到标准目录格式： - 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>
This commit is contained in:
@@ -0,0 +1,938 @@
|
||||
"""
|
||||
概率与不确定性示例 (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()
|
||||
Reference in New Issue
Block a user