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>
1043 lines
32 KiB
Python
1043 lines
32 KiB
Python
"""
|
|
不确定性分析示例 (Uncertainty Analysis Example)
|
|
==============================================
|
|
|
|
本示例展示空间智能系统中的不确定性分析方法。
|
|
在空间决策中，不确定性来自数据、模型和参数等多个方面。
|
|
|
|
核心概念:
|
|
1. 不确定性来源 - 数据误差、模型简化、参数变异
|
|
2. 不确定性传播 - 输入不确定性如何影响输出
|
|
3. 蒙特卡洛分析 - 随机采样评估不确定性
|
|
4. 敏感性分析 - 识别关键不确定性源
|
|
5. 场景分析 - 不同假设下的结果比较
|
|
|
|
应用场景:
|
|
- 风险评估
|
|
- 决策稳健性分析
|
|
- 模型可信度评估
|
|
- 数据质量评估
|
|
|
|
作者: CC4SI 项目组
|
|
"""
|
|
|
|
import math
|
|
import random
|
|
from typing import List, Dict, Tuple, Optional, Any, Callable
|
|
from dataclasses import dataclass, field
|
|
from enum import Enum
|
|
from abc import ABC, abstractmethod
|
|
import statistics
|
|
|
|
|
|
# ============================================================================
|
|
# 不确定性类型
|
|
# ============================================================================
|
|
|
|
class UncertaintyType(Enum):
|
|
"""不确定性类型"""
|
|
EPISTEMIC = "epistemic" # 认识不确定性 (可通过更多知识减少)
|
|
ALEATORY = "aleatory" # 偶然不确定性 (固有随机性)
|
|
PARAMETRIC = "parametric" # 参数不确定性
|
|
STRUCTURAL = "structural" # 结构不确定性 (模型形式)
|
|
DATA = "data" # 数据不确定性
|
|
|
|
|
|
@dataclass
|
|
class UncertainValue:
|
|
"""
|
|
不确定值
|
|
|
|
表示一个带有不确定性的数值。
|
|
"""
|
|
value: float
|
|
uncertainty: float # 标准差或误差范围
|
|
uncertainty_type: UncertaintyType = UncertaintyType.EPISTEMIC
|
|
distribution: str = "normal" # 假设的分布类型
|
|
|
|
@property
|
|
def coefficient_of_variation(self) -> float:
|
|
"""变异系数"""
|
|
if self.value == 0:
|
|
return float('inf')
|
|
return self.uncertainty / abs(self.value)
|
|
|
|
def confidence_interval(self, confidence: float = 0.95) -> Tuple[float, float]:
|
|
"""
|
|
计算置信区间
|
|
|
|
Args:
|
|
confidence: 置信水平
|
|
|
|
Returns:
|
|
(下界, 上界)
|
|
"""
|
|
if self.distribution == "normal":
|
|
# 使用正态分布
|
|
z_scores = {0.90: 1.645, 0.95: 1.96, 0.99: 2.576}
|
|
z = z_scores.get(confidence, 1.96)
|
|
return (
|
|
self.value - z * self.uncertainty,
|
|
self.value + z * self.uncertainty
|
|
)
|
|
else:
|
|
# 简单区间
|
|
return (
|
|
self.value - self.uncertainty,
|
|
self.value + self.uncertainty
|
|
)
|
|
|
|
def sample(self) -> float:
|
|
"""从分布中采样"""
|
|
if self.distribution == "normal":
|
|
return random.gauss(self.value, self.uncertainty)
|
|
elif self.distribution == "uniform":
|
|
return random.uniform(
|
|
self.value - self.uncertainty,
|
|
self.value + self.uncertainty
|
|
)
|
|
else:
|
|
return self.value
|
|
|
|
def __repr__(self) -> str:
|
|
return f"{self.value:.2f} ± {self.uncertainty:.2f}"
|
|
|
|
|
|
# ============================================================================
|
|
# 概率分布
|
|
# ============================================================================
|
|
|
|
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
|
|
|
|
|
|
class NormalDistribution(ProbabilityDistribution):
|
|
"""正态分布"""
|
|
|
|
def __init__(self, mu: float, sigma: float, name: str = ""):
|
|
super().__init__(name)
|
|
self.mu = mu
|
|
self.sigma = sigma
|
|
|
|
def sample(self) -> float:
|
|
return random.gauss(self.mu, self.sigma)
|
|
|
|
def mean(self) -> float:
|
|
return self.mu
|
|
|
|
def std(self) -> float:
|
|
return self.sigma
|
|
|
|
|
|
class UniformDistribution(ProbabilityDistribution):
|
|
"""均匀分布"""
|
|
|
|
def __init__(self, a: float, b: float, name: str = ""):
|
|
super().__init__(name)
|
|
self.a = a
|
|
self.b = b
|
|
|
|
def sample(self) -> float:
|
|
return random.uniform(self.a, self.b)
|
|
|
|
def mean(self) -> float:
|
|
return (self.a + self.b) / 2
|
|
|
|
def std(self) -> float:
|
|
return (self.b - self.a) / math.sqrt(12)
|
|
|
|
|
|
class TriangularDistribution(ProbabilityDistribution):
|
|
"""三角分布"""
|
|
|
|
def __init__(self, a: float, b: float, c: float, name: str = ""):
|
|
"""
|
|
三角分布
|
|
|
|
Args:
|
|
a: 最小值
|
|
b: 最大值
|
|
c: 众数
|
|
"""
|
|
super().__init__(name)
|
|
self.a = a
|
|
self.b = b
|
|
self.c = c
|
|
|
|
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)
|
|
|
|
|
|
# ============================================================================
|
|
# 不确定性传播
|
|
# ============================================================================
|
|
|
|
class UncertaintyPropagator:
|
|
"""
|
|
不确定性传播器
|
|
|
|
分析输入不确定性如何影响输出。
|
|
"""
|
|
|
|
def __init__(self, model: Callable[[Dict[str, float]], float]):
|
|
"""
|
|
初始化传播器
|
|
|
|
Args:
|
|
model: 输入字典到输出值的函数
|
|
"""
|
|
self.model = model
|
|
|
|
def first_order_second_moment(self,
|
|
inputs: Dict[str, UncertainValue]) -> UncertainValue:
|
|
"""
|
|
一阶二矩法 (FOSM)
|
|
|
|
使用一阶泰勒展开近似传播不确定性。
|
|
|
|
Args:
|
|
inputs: 不确定输入字典
|
|
|
|
Returns:
|
|
不确定输出
|
|
"""
|
|
# 计算名义值
|
|
nominal_inputs = {k: v.value for k, v in inputs.items()}
|
|
nominal_output = self.model(nominal_inputs)
|
|
|
|
# 计算灵敏度 (数值微分)
|
|
sensitivities = {}
|
|
epsilon = 1e-6
|
|
|
|
for name, uncertain_val in inputs.items():
|
|
perturbed = nominal_inputs.copy()
|
|
perturbed[name] += epsilon
|
|
output_plus = self.model(perturbed)
|
|
sensitivity = (output_plus - nominal_output) / epsilon
|
|
sensitivities[name] = sensitivity
|
|
|
|
# 计算输出方差 (假设输入独立)
|
|
output_variance = 0.0
|
|
for name, uncertain_val in inputs.items():
|
|
sensitivity = sensitivities[name]
|
|
output_variance += (sensitivity * uncertain_val.uncertainty) ** 2
|
|
|
|
output_std = math.sqrt(output_variance)
|
|
|
|
return UncertainValue(
|
|
value=nominal_output,
|
|
uncertainty=output_std,
|
|
uncertainty_type=UncertaintyType.EPISTEMIC
|
|
)
|
|
|
|
def monte_carlo_propagation(self,
|
|
input_distributions: Dict[str, ProbabilityDistribution],
|
|
n_samples: int = 10000) -> Dict[str, Any]:
|
|
"""
|
|
蒙特卡洛传播
|
|
|
|
Args:
|
|
input_distributions: 输入分布字典
|
|
n_samples: 采样次数
|
|
|
|
Returns:
|
|
统计结果字典
|
|
"""
|
|
samples = []
|
|
|
|
for _ in range(n_samples):
|
|
# 采样输入
|
|
inputs = {name: dist.sample() for name, dist in input_distributions.items()}
|
|
# 计算输出
|
|
output = self.model(inputs)
|
|
samples.append(output)
|
|
|
|
# 计算统计量
|
|
return {
|
|
"mean": statistics.mean(samples),
|
|
"std": statistics.stdev(samples) if len(samples) > 1 else 0,
|
|
"min": min(samples),
|
|
"max": max(samples),
|
|
"median": statistics.median(samples),
|
|
"percentiles": {
|
|
5: self._percentile(samples, 5),
|
|
25: self._percentile(samples, 25),
|
|
75: self._percentile(samples, 75),
|
|
95: self._percentile(samples, 95)
|
|
},
|
|
"samples": samples
|
|
}
|
|
|
|
def _percentile(self, data: List[float], p: float) -> float:
|
|
"""计算百分位数"""
|
|
sorted_data = sorted(data)
|
|
index = int(p / 100 * len(sorted_data))
|
|
return sorted_data[min(index, len(sorted_data) - 1)]
|
|
|
|
|
|
# ============================================================================
|
|
# 敏感性分析
|
|
# ============================================================================
|
|
|
|
@dataclass
|
|
class SensitivityMeasure:
|
|
"""敏感性度量"""
|
|
parameter_name: str
|
|
sensitivity: float # 敏感性系数
|
|
rank: int = 0 # 排名
|
|
method: str = "" # 计算方法
|
|
|
|
|
|
class SensitivityAnalyzer:
|
|
"""
|
|
敏感性分析器
|
|
|
|
识别对输出影响最大的输入参数。
|
|
"""
|
|
|
|
def __init__(self, model: Callable[[Dict[str, float]], float]):
|
|
"""
|
|
初始化分析器
|
|
|
|
Args:
|
|
model: 输入字典到输出值的函数
|
|
"""
|
|
self.model = model
|
|
|
|
def local_sensitivity(self,
|
|
nominal_values: Dict[str, float],
|
|
perturbation: float = 0.01) -> List[SensitivityMeasure]:
|
|
"""
|
|
局部敏感性分析
|
|
|
|
Args:
|
|
nominal_values: 名义输入值
|
|
perturbation: 扰动比例
|
|
|
|
Returns:
|
|
敏感性度量列表
|
|
"""
|
|
nominal_output = self.model(nominal_values)
|
|
sensitivities = []
|
|
|
|
for param_name in nominal_values.keys():
|
|
# 正向扰动
|
|
perturbed_plus = nominal_values.copy()
|
|
perturbed_plus[param_name] *= (1 + perturbation)
|
|
output_plus = self.model(perturbed_plus)
|
|
|
|
# 计算敏感性 (归一化)
|
|
delta_output = output_plus - nominal_output
|
|
sensitivity = delta_output / (perturbation * nominal_values[param_name])
|
|
sensitivities.append({
|
|
"parameter": param_name,
|
|
"sensitivity": sensitivity,
|
|
"delta_output": delta_output
|
|
})
|
|
|
|
# 排序 (绝对值)
|
|
sensitivities.sort(key=lambda x: abs(x["sensitivity"]), reverse=True)
|
|
|
|
measures = []
|
|
for i, s in enumerate(sensitivities, 1):
|
|
measures.append(SensitivityMeasure(
|
|
parameter_name=s["parameter"],
|
|
sensitivity=s["sensitivity"],
|
|
rank=i,
|
|
method="local"
|
|
))
|
|
|
|
return measures
|
|
|
|
def variance_based_sensitivity(self,
|
|
input_distributions: Dict[str, ProbabilityDistribution],
|
|
n_samples: int = 1000) -> List[SensitivityMeasure]:
|
|
"""
|
|
基于方差的敏感性分析 (Sobol指数近似)
|
|
|
|
Args:
|
|
input_distributions: 输入分布
|
|
n_samples: 样本数
|
|
|
|
Returns:
|
|
敏感性度量列表
|
|
"""
|
|
# 使用随机平衡设计近似一阶Sobol指数
|
|
samples = []
|
|
for _ in range(n_samples):
|
|
inputs = {name: dist.sample() for name, dist in input_distributions.items()}
|
|
outputs = self.model(inputs)
|
|
samples.append((inputs, outputs))
|
|
|
|
# 计算总方差
|
|
output_values = [s[1] for s in samples]
|
|
total_variance = statistics.variance(output_values) if len(output_values) > 1 else 0
|
|
|
|
if total_variance == 0:
|
|
return []
|
|
|
|
sensitivities = []
|
|
|
|
for param_name in input_distributions.keys():
|
|
# 计算条件方差 (简化: 使用回归方法)
|
|
param_values = [s[0][param_name] for s in samples]
|
|
|
|
# 线性回归 R²
|
|
mean_x = statistics.mean(param_values)
|
|
mean_y = statistics.mean(output_values)
|
|
|
|
numerator = sum((x - mean_x) * (y - mean_y)
|
|
for x, y in zip(param_values, output_values))
|
|
denominator_x = sum((x - mean_x)**2 for x in param_values)
|
|
denominator_y = sum((y - mean_y)**2 for y in output_values)
|
|
|
|
if denominator_x == 0 or denominator_y == 0:
|
|
first_order = 0
|
|
else:
|
|
correlation = numerator / math.sqrt(denominator_x * denominator_y)
|
|
first_order = correlation ** 2
|
|
|
|
sensitivities.append({
|
|
"parameter": param_name,
|
|
"sensitivity": first_order,
|
|
"variance_contribution": first_order * total_variance
|
|
})
|
|
|
|
# 排序
|
|
sensitivities.sort(key=lambda x: x["sensitivity"], reverse=True)
|
|
|
|
measures = []
|
|
for i, s in enumerate(sensitivities, 1):
|
|
measures.append(SensitivityMeasure(
|
|
parameter_name=s["parameter"],
|
|
sensitivity=s["sensitivity"],
|
|
rank=i,
|
|
method="sobol"
|
|
))
|
|
|
|
return measures
|
|
|
|
|
|
# ============================================================================
|
|
# 场景分析
|
|
# ============================================================================
|
|
|
|
@dataclass
|
|
class Scenario:
|
|
"""场景"""
|
|
name: str
|
|
description: str
|
|
parameters: Dict[str, float]
|
|
probability: float = 1.0
|
|
|
|
|
|
class ScenarioAnalyzer:
|
|
"""
|
|
场景分析器
|
|
|
|
评估不同假设情景下的结果。
|
|
"""
|
|
|
|
def __init__(self, model: Callable[[Dict[str, float]], float]):
|
|
"""
|
|
初始化分析器
|
|
|
|
Args:
|
|
model: 输入字典到输出值的函数
|
|
"""
|
|
self.model = model
|
|
|
|
def evaluate_scenarios(self,
|
|
scenarios: List[Scenario]) -> List[Dict[str, Any]]:
|
|
"""
|
|
评估多个场景
|
|
|
|
Args:
|
|
scenarios: 场景列表
|
|
|
|
Returns:
|
|
评估结果列表
|
|
"""
|
|
results = []
|
|
|
|
for scenario in scenarios:
|
|
output = self.model(scenario.parameters)
|
|
|
|
results.append({
|
|
"scenario": scenario.name,
|
|
"description": scenario.description,
|
|
"probability": scenario.probability,
|
|
"output": output,
|
|
"parameters": scenario.parameters.copy()
|
|
})
|
|
|
|
return results
|
|
|
|
def generate_scenarios(self,
|
|
base_parameters: Dict[str, float],
|
|
variations: Dict[str, Tuple[float, float]]) -> List[Scenario]:
|
|
"""
|
|
生成场景 (基准、乐观、悲观)
|
|
|
|
Args:
|
|
base_parameters: 基准参数
|
|
variations: 参数变化范围字典
|
|
|
|
Returns:
|
|
场景列表
|
|
"""
|
|
scenarios = [
|
|
Scenario(
|
|
name="基准",
|
|
description="预期情况",
|
|
parameters=base_parameters.copy(),
|
|
probability=0.5
|
|
)
|
|
]
|
|
|
|
# 乐观场景
|
|
optimistic = base_parameters.copy()
|
|
for param, (low, high) in variations.items():
|
|
if param in base_parameters:
|
|
# 选择有利方向的值
|
|
if "成本" in param or "cost" in param.lower():
|
|
optimistic[param] = low # 成本取低
|
|
else:
|
|
optimistic[param] = high # 其他取高
|
|
|
|
scenarios.append(Scenario(
|
|
name="乐观",
|
|
description="最佳情况",
|
|
parameters=optimistic,
|
|
probability=0.25
|
|
))
|
|
|
|
# 悲观场景
|
|
pessimistic = base_parameters.copy()
|
|
for param, (low, high) in variations.items():
|
|
if param in base_parameters:
|
|
if "成本" in param or "cost" in param.lower():
|
|
pessimistic[param] = high # 成本取高
|
|
else:
|
|
pessimistic[param] = low # 其他取低
|
|
|
|
scenarios.append(Scenario(
|
|
name="悲观",
|
|
description="最差情况",
|
|
parameters=pessimistic,
|
|
probability=0.25
|
|
))
|
|
|
|
return scenarios
|
|
|
|
|
|
# ============================================================================
|
|
# 稳健性分析
|
|
# ============================================================================
|
|
|
|
class RobustnessAnalyzer:
|
|
"""
|
|
稳健性分析器
|
|
|
|
评估决策在不同条件下的稳健程度。
|
|
"""
|
|
|
|
def __init__(self, model: Callable[[Dict[str, float]], float]):
|
|
"""
|
|
初始化分析器
|
|
|
|
Args:
|
|
model: 输入字典到输出值的函数
|
|
"""
|
|
self.model = model
|
|
|
|
def worst_case_analysis(self,
|
|
nominal_values: Dict[str, float],
|
|
uncertainties: Dict[str, float],
|
|
n_samples: int = 1000) -> Dict[str, Any]:
|
|
"""
|
|
最坏情况分析
|
|
|
|
Args:
|
|
nominal_values: 名义值
|
|
uncertainties: 不确定性范围
|
|
n_samples: 采样次数
|
|
|
|
Returns:
|
|
分析结果
|
|
"""
|
|
samples = []
|
|
|
|
for _ in range(n_samples):
|
|
perturbed = {}
|
|
for param, nominal in nominal_values.items():
|
|
uncertainty = uncertainties.get(param, 0)
|
|
# 均匀采样
|
|
perturbed[param] = random.uniform(
|
|
nominal - uncertainty,
|
|
nominal + uncertainty
|
|
)
|
|
|
|
output = self.model(perturbed)
|
|
samples.append(output)
|
|
|
|
return {
|
|
"nominal": self.model(nominal_values),
|
|
"best": max(samples),
|
|
"worst": min(samples),
|
|
"range": max(samples) - min(samples),
|
|
"mean": statistics.mean(samples),
|
|
"std": statistics.stdev(samples) if len(samples) > 1 else 0,
|
|
"percentile_5": self._percentile(samples, 5),
|
|
"percentile_95": self._percentile(samples, 95)
|
|
}
|
|
|
|
def regret_analysis(self,
|
|
alternatives: List[Dict[str, float]],
|
|
scenarios: List[Dict[str, float]],
|
|
alternative_names: List[str] = None) -> Dict[str, Any]:
|
|
"""
|
|
后悔值分析
|
|
|
|
Args:
|
|
alternatives: 备选方案列表
|
|
scenarios: 场景列表
|
|
alternative_names: 方案名称
|
|
|
|
Returns:
|
|
后悔值分析结果
|
|
"""
|
|
if alternative_names is None:
|
|
alternative_names = [f"方案{i+1}" for i in range(len(alternatives))]
|
|
|
|
# 计算每个方案在各场景下的结果
|
|
results_matrix = []
|
|
for alt in alternatives:
|
|
row = []
|
|
for scenario in scenarios:
|
|
# 合并参数
|
|
combined = {**alt, **scenario}
|
|
output = self.model(combined)
|
|
row.append(output)
|
|
results_matrix.append(row)
|
|
|
|
# 找出每个场景下的最优结果
|
|
best_per_scenario = []
|
|
for j in range(len(scenarios)):
|
|
column = [results_matrix[i][j] for i in range(len(alternatives))]
|
|
best_per_scenario.append(max(column) if column else 0)
|
|
|
|
# 计算后悔值矩阵
|
|
regret_matrix = []
|
|
for i in range(len(alternatives)):
|
|
regrets = []
|
|
for j in range(len(scenarios)):
|
|
regret = best_per_scenario[j] - results_matrix[i][j]
|
|
regrets.append(regret)
|
|
regret_matrix.append(regrets)
|
|
|
|
# 计算最大后悔值
|
|
max_regrets = [max(regrets) for regrets in regret_matrix]
|
|
|
|
# 排序
|
|
sorted_indices = sorted(range(len(max_regrets)), key=lambda i: max_regrets[i])
|
|
|
|
return {
|
|
"results_matrix": results_matrix,
|
|
"regret_matrix": regret_matrix,
|
|
"max_regrets": max_regrets,
|
|
"minimax_regret_choice": sorted_indices[0] if sorted_indices else None,
|
|
"ranking": sorted_indices,
|
|
"alternative_names": alternative_names
|
|
}
|
|
|
|
def _percentile(self, data: List[float], p: float) -> float:
|
|
"""计算百分位数"""
|
|
sorted_data = sorted(data)
|
|
index = int(p / 100 * len(sorted_data))
|
|
return sorted_data[min(index, len(sorted_data) - 1)]
|
|
|
|
|
|
# ============================================================================
|
|
# 空间不确定性应用
|
|
# ============================================================================
|
|
|
|
class SpatialUncertaintyModel:
|
|
"""
|
|
空间不确定性模型
|
|
|
|
处理空间数据中的不确定性。
|
|
"""
|
|
|
|
@staticmethod
|
|
def uncertain_distance(p1: Tuple[float, float],
|
|
p2: Tuple[float, float],
|
|
position_error: float = 5.0) -> float:
|
|
"""
|
|
带不确定性的距离计算
|
|
|
|
Args:
|
|
p1: 点1坐标
|
|
p2: 点2坐标
|
|
position_error: 位置误差标准差
|
|
|
|
Returns:
|
|
采样距离
|
|
"""
|
|
# 添加位置误差
|
|
x1 = p1[0] + random.gauss(0, position_error)
|
|
y1 = p1[1] + random.gauss(0, position_error)
|
|
x2 = p2[0] + random.gauss(0, position_error)
|
|
y2 = p2[1] + random.gauss(0, position_error)
|
|
|
|
return math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
|
|
|
|
@staticmethod
|
|
def uncertain_interpolation(target_point: Tuple[float, float],
|
|
sample_points: List[Tuple[Tuple[float, float], float]],
|
|
measurement_error: float = 0.1,
|
|
power: float = 2.0) -> float:
|
|
"""
|
|
带不确定性的空间插值 (IDW)
|
|
|
|
Args:
|
|
target_point: 目标点
|
|
sample_points: 样本点列表 [(坐标, 值), ...]
|
|
measurement_error: 测量误差标准差
|
|
power: IDW幂次
|
|
|
|
Returns:
|
|
插值结果
|
|
"""
|
|
tx, ty = target_point
|
|
|
|
numerator = 0.0
|
|
denominator = 0.0
|
|
|
|
for (sx, sy), value in sample_points:
|
|
# 计算距离
|
|
distance = math.sqrt((tx - sx)**2 + (ty - sy)**2)
|
|
|
|
if distance < 1e-10:
|
|
# 几乎重合，返回该值加误差
|
|
return value + random.gauss(0, measurement_error)
|
|
|
|
# 添加测量误差
|
|
observed_value = value + random.gauss(0, measurement_error)
|
|
|
|
# 计算权重
|
|
weight = 1.0 / (distance ** power)
|
|
|
|
numerator += weight * observed_value
|
|
denominator += weight
|
|
|
|
if denominator == 0:
|
|
return 0.0
|
|
|
|
return numerator / denominator
|
|
|
|
|
|
# ============================================================================
|
|
# 主程序
|
|
# ========================================================================
|
|
|
|
def main():
|
|
"""主程序 - 演示不确定性分析的使用"""
|
|
|
|
print("="*70)
|
|
print("不确定性分析示例演示")
|
|
print("="*70)
|
|
|
|
random.seed(42)
|
|
|
|
# ========================================================================
|
|
# 1. 不确定值表示
|
|
# ========================================================================
|
|
print("\n[部分 1] 不确定值表示")
|
|
print("-" * 50)
|
|
|
|
# 创建不确定值
|
|
population = UncertainValue(
|
|
value=10000,
|
|
uncertainty=500,
|
|
uncertainty_type=UncertaintyType.EPISTEMIC,
|
|
distribution="normal"
|
|
)
|
|
|
|
print(f"\n人口估计: {population}")
|
|
print(f" 变异系数: {population.coefficient_of_variation:.3f}")
|
|
print(f" 95% 置信区间: {population.confidence_interval(0.95)}")
|
|
|
|
# 采样演示
|
|
print(f"\n采样示例:")
|
|
samples = [population.sample() for _ in range(5)]
|
|
print(f" {samples}")
|
|
|
|
# ========================================================================
|
|
# 2. 不确定性传播
|
|
# ========================================================================
|
|
print("\n\n[部分 2] 不确定性传播")
|
|
print("-" * 50)
|
|
|
|
# 定义模型: 地块价值 = 面积 * 单价 - 开发成本
|
|
def land_value_model(inputs):
|
|
area = inputs["area"]
|
|
unit_price = inputs["unit_price"]
|
|
development_cost = inputs["development_cost"]
|
|
return area * unit_price - development_cost
|
|
|
|
propagator = UncertaintyPropagator(land_value_model)
|
|
|
|
# 定义不确定输入
|
|
uncertain_inputs = {
|
|
"area": UncertainValue(1000, 50), # 面积: 1000 ± 50 平方米
|
|
"unit_price": UncertainValue(5000, 300), # 单价: 5000 ± 300 元/平方米
|
|
"development_cost": UncertainValue(100000, 10000) # 开发成本
|
|
}
|
|
|
|
print("\n输入不确定性:")
|
|
for name, val in uncertain_inputs.items():
|
|
print(f" {name}: {val}")
|
|
|
|
# 一阶二矩法
|
|
fosm_result = propagator.first_order_second_moment(uncertain_inputs)
|
|
print(f"\n一阶二矩法 (FOSM) 结果:")
|
|
print(f" 地块价值: {fosm_result}")
|
|
print(f" 95% 置信区间: {fosm_result.confidence_interval(0.95)}")
|
|
|
|
# 蒙特卡洛传播
|
|
input_dists = {
|
|
"area": NormalDistribution(1000, 50),
|
|
"unit_price": NormalDistribution(5000, 300),
|
|
"development_cost": NormalDistribution(100000, 10000)
|
|
}
|
|
|
|
mc_result = propagator.monte_carlo_propagation(input_dists, n_samples=10000)
|
|
print(f"\n蒙特卡洛传播结果:")
|
|
print(f" 均值: {mc_result['mean']:.0f}")
|
|
print(f" 标准差: {mc_result['std']:.0f}")
|
|
print(f" 范围: [{mc_result['min']:.0f}, {mc_result['max']:.0f}]")
|
|
print(f" 90% 置信区间: [{mc_result['percentiles'][5]:.0f}, "
|
|
f"{mc_result['percentiles'][95]:.0f}]")
|
|
|
|
# ========================================================================
|
|
# 3. 敏感性分析
|
|
# ========================================================================
|
|
print("\n\n[部分 3] 敏感性分析")
|
|
print("-" * 50)
|
|
|
|
# 局部敏感性
|
|
analyzer = SensitivityAnalyzer(land_value_model)
|
|
|
|
nominal_values = {
|
|
"area": 1000,
|
|
"unit_price": 5000,
|
|
"development_cost": 100000
|
|
}
|
|
|
|
local_sens = analyzer.local_sensitivity(nominal_values, perturbation=0.01)
|
|
|
|
print("\n局部敏感性分析:")
|
|
print(f"{'排名':<6} {'参数':<15} {'敏感性系数':<15} {'影响':<15}")
|
|
print("-" * 55)
|
|
for s in local_sens:
|
|
impact = "高" if abs(s.sensitivity) > 100 else "中" if abs(s.sensitivity) > 10 else "低"
|
|
print(f"{s.rank:<6} {s.parameter_name:<15} {s.sensitivity:<15.2f} {impact:<15}")
|
|
|
|
# 方差基敏感性
|
|
var_sens = analyzer.variance_based_sensitivity(input_dists, n_samples=1000)
|
|
|
|
print("\n基于方差的敏感性 (Sobol指数近似):")
|
|
print(f"{'排名':<6} {'参数':<15} {'一阶效应':<15} {'贡献':<15}")
|
|
print("-" * 55)
|
|
for s in var_sens:
|
|
contrib = f"{s.sensitivity * 100:.1f}%"
|
|
print(f"{s.rank:<6} {s.parameter_name:<15} {s.sensitivity:<15.4f} {contrib:<15}")
|
|
|
|
# ========================================================================
|
|
# 4. 场景分析
|
|
# ========================================================================
|
|
print("\n\n[部分 4] 场景分析 - 房地产项目评估")
|
|
print("-" * 50)
|
|
|
|
# 定义项目评估模型
|
|
def project_evaluation(inputs):
|
|
revenue = inputs["area"] * inputs["selling_price"]
|
|
cost = inputs["land_cost"] + inputs["construction_cost"] * inputs["area"]
|
|
return revenue - cost
|
|
|
|
scenario_analyzer = ScenarioAnalyzer(project_evaluation)
|
|
|
|
base_params = {
|
|
"area": 10000,
|
|
"selling_price": 15000,
|
|
"land_cost": 50000000,
|
|
"construction_cost": 8000
|
|
}
|
|
|
|
variations = {
|
|
"area": (8000, 12000),
|
|
"selling_price": (12000, 18000),
|
|
"land_cost": (40000000, 60000000),
|
|
"construction_cost": (7000, 9000)
|
|
}
|
|
|
|
scenarios = scenario_analyzer.generate_scenarios(base_params, variations)
|
|
|
|
print("\n生成的场景:")
|
|
for s in scenarios:
|
|
print(f" {s.name}: {s.description} (概率: {s.probability})")
|
|
|
|
# 评估场景
|
|
results = scenario_analyzer.evaluate_scenarios(scenarios)
|
|
|
|
print("\n场景评估结果:")
|
|
print(f"{'场景':<10} {'利润(万元)':<15} {'概率':<10}")
|
|
print("-" * 40)
|
|
for r in results:
|
|
profit_wan = r["output"] / 10000
|
|
print(f"{r['scenario']:<10} {profit_wan:<15.1f} {r['probability']:<10.1%}")
|
|
|
|
# 期望值
|
|
expected_value = sum(r["output"] * r["probability"] for r in results)
|
|
print(f"\n期望利润: {expected_value / 10000:.1f} 万元")
|
|
|
|
# ========================================================================
|
|
# 5. 稳健性分析
|
|
# ========================================================================
|
|
print("\n\n[部分 5] 稳健性分析")
|
|
print("-" * 50)
|
|
|
|
robustness = RobustnessAnalyzer(project_evaluation)
|
|
|
|
# 最坏情况分析
|
|
uncertainties = {
|
|
"area": 1000,
|
|
"selling_price": 2000,
|
|
"land_cost": 5000000,
|
|
"construction_cost": 500
|
|
}
|
|
|
|
worst_case = robustness.worst_case_analysis(base_params, uncertainties, n_samples=1000)
|
|
|
|
print("\n最坏情况分析:")
|
|
print(f" 名义利润: {worst_case['nominal'] / 10000:.1f} 万元")
|
|
print(f" 最好情况: {worst_case['best'] / 10000:.1f} 万元")
|
|
print(f" 最差情况: {worst_case['worst'] / 10000:.1f} 万元")
|
|
print(f" 变化范围: {worst_case['range'] / 10000:.1f} 万元")
|
|
print(f" 90% 置信区间: [{worst_case['percentile_5'] / 10000:.1f}, "
|
|
f"{worst_case['percentile_95'] / 10000:.1f}] 万元")
|
|
|
|
# 后悔值分析
|
|
alternatives = [
|
|
{"area": 8000, "selling_price": 15000, "land_cost": 40000000, "construction_cost": 8000},
|
|
{"area": 10000, "selling_price": 15000, "land_cost": 50000000, "construction_cost": 8000},
|
|
{"area": 12000, "selling_price": 15000, "land_cost": 60000000, "construction_cost": 8000},
|
|
]
|
|
|
|
scenarios_list = [
|
|
{"selling_price": 13000}, # 价格下跌
|
|
{"selling_price": 15000}, # 价格平稳
|
|
{"selling_price": 17000}, # 价格上涨
|
|
]
|
|
|
|
regret_result = robustness.regret_analysis(
|
|
alternatives, scenarios_list,
|
|
alternative_names=["小规模", "中规模", "大规模"]
|
|
)
|
|
|
|
print("\n后悔值分析:")
|
|
print(f"{'方案':<10} {'最大后悔值(万元)':<20}")
|
|
print("-" * 35)
|
|
for i, name in enumerate(regret_result["alternative_names"]):
|
|
max_regret_wan = regret_result["max_regrets"][i] / 10000
|
|
print(f"{name:<10} {max_regret_wan:<20.1f}")
|
|
|
|
minimax_choice = regret_result["minimax_regret_choice"]
|
|
if minimax_choice is not None:
|
|
best_name = regret_result["alternative_names"][minimax_choice]
|
|
print(f"\n最小最大后悔值推荐: {best_name}")
|
|
|
|
# ========================================================================
|
|
# 6. 空间不确定性
|
|
# ========================================================================
|
|
print("\n\n[部分 6] 空间不确定性应用")
|
|
print("-" * 50)
|
|
|
|
# GPS定位不确定下的距离测量
|
|
point_a = (100, 200)
|
|
point_b = (150, 250)
|
|
|
|
print("\n带位置误差的距离测量:")
|
|
print(f" 点A: {point_a}")
|
|
print(f" 点B: {point_b}")
|
|
print(f" 理论距离: {math.sqrt((150-100)**2 + (250-200)**2):.2f}")
|
|
|
|
# 多次测量
|
|
measurements = [SpatialUncertaintyModel.uncertain_distance(
|
|
point_a, point_b, position_error=5
|
|
) for _ in range(10)]
|
|
|
|
print(f" 实际测量 (10次): {[f'{m:.1f}' for m in measurements]}")
|
|
print(f" 平均: {statistics.mean(measurements):.2f} ± {statistics.stdev(measurements):.2f}")
|
|
|
|
# 不确定插值
|
|
print("\n带测量误差的空间插值:")
|
|
samples = [
|
|
((0, 0), 10),
|
|
((100, 0), 20),
|
|
((0, 100), 15),
|
|
((100, 100), 25)
|
|
]
|
|
target = (50, 50)
|
|
|
|
interpolated_values = [
|
|
SpatialUncertaintyModel.uncertain_interpolation(
|
|
target, samples, measurement_error=1
|
|
) for _ in range(10)
|
|
]
|
|
|
|
print(f" 目标点: {target}")
|
|
print(f" 插值结果 (10次): {[f'{v:.1f}' for v in interpolated_values]}")
|
|
print(f" 平均: {statistics.mean(interpolated_values):.2f} ± "
|
|
f"{statistics.stdev(interpolated_values):.2f}")
|
|
|
|
print("\n" + "="*70)
|
|
print("演示完成!")
|
|
print("="*70)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|