使用MLE拟合对数正态分布构建易损性曲线的优化问题

问题概述 #

针对不同峰值地面加速度（PGA）下的损伤概率数据，使用最大似然估计（MLE）拟合对数正态分布时出现以下问题：

• 轻微损伤等级：优化未收敛，参数估计无实际意义
• 严重损伤等级：拟合成功，结果合理
• 倒塌损伤等级：达到最大迭代次数，参数估计异常

问题根源分析 #

1. 似然函数数值不稳定：直接计算二项分布PMF后取对数，当PMF为0时会得到-inf，代码中错误添加1e+10导致似然和异常，最终出现inf值。
2. 初始参数选择不合理：所有损伤等级都使用(1,1)作为初始值，与轻微/倒塌损伤的物理规律（轻微损伤低PGA概率高、倒塌损伤高PGA概率高）偏差过大，无梯度优化算法易陷入局部最优。
3. 参数约束缺失：未限制sigma（对数正态中位值）和beta（对数标准差）为正数，导致倒塌损伤的sigma出现极大值。
4. 冗余打印干扰：迭代过程中频繁打印参数，减慢优化速度并可能影响数值稳定性。

修正方案 #

1. 修复似然函数数值稳定性 #

使用scipy.stats.binom.logpmf直接计算对数似然，避免先计算PMF再取对数的数值问题，移除错误的偏移量操作。

2. 添加参数约束 #

采用支持边界约束的L-BFGS-B优化方法，设置sigma和beta的下界为1e-6，确保参数为正数。

3. 优化初始参数 #

根据损伤等级的物理规律设置初始值：

• 轻微损伤：sigma=0.2（低中位PGA），beta=0.5
• 严重损伤：sigma=0.4，beta=0.5
• 倒塌损伤：sigma=0.5（高中位PGA），beta=0.5

4. 移除冗余打印 #

删除似然函数内的打印语句，减少优化干扰。

修正后的完整代码 #

from functools import partial
import numpy as np
from scipy import optimize, stats
import matplotlib.pyplot as plt

PGA_Values = np.array([0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40,
               0.42, 0.44, 0.46, 0.48, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, 0.64, 0.66, 0.68, 0.70,
               0.72])

# 轻微损伤数据
slight_damage_analyses = np.array([46, 11, 18, 54, 30, 1482, 17, 73, 1082, 341, 798, 9509, 4226, 107, 7, 1,1,
                               3, 11, 103, 9, 27, 21, 22, 40, 22, 90, 98, 1006, 402, 1]) 
slight_damage_collapses = np.array([29, 5, 12, 33, 18, 1314, 10, 47, 656, 217, 611, 6751, 2734, 69, 2, 0, 0,
                                0, 5, 49, 5, 0, 7, 8, 13, 3, 55, 67, 342, 260, 1])

# 严重损伤数据
heavily_damage_analyses = np.array([46, 11, 18, 54, 30, 1482, 17, 73, 1082, 341, 798, 9509, 4226, 107, 7, 1,1,
                               3, 11, 103, 9, 27, 21, 22, 40, 22, 90, 98, 1006, 402, 1]) 
heavily_damage_collapses = np.array([15, 6, 3, 18, 12, 150, 5, 18, 342, 83, 160, 2244, 1107, 27, 1, 1, 1, 1,
                                4, 49, 3, 6, 10, 12, 17, 6, 34, 27, 597, 114, 0])

# 倒塌损伤数据
collapsed_damage_analyses = np.array([46, 11, 18, 54, 30, 1482, 17, 73, 1082, 341, 798, 9509, 4226, 107, 7, 1,1,
                               3, 11, 103, 9, 27, 21, 22, 40, 22, 90, 98, 1006, 402, 1]) 
collapsed_damage_collapses = np.array([2, 0, 3, 3, 0, 18, 2, 8, 84, 41, 27, 514, 385, 11, 4, 1, 0, 2, 2, 5, 1,
                                21, 4, 2, 10, 13, 1, 4, 67, 28, 0])

def neg_log_likelihood_sum(params, im_l, no_a, no_c):
    sigma, beta = params
    # 对数正态CDF等价于正态分布在log(PGA)处的CDF，mu=log(sigma)
    theoretical_prob = stats.norm(loc=np.log(sigma), scale=beta).cdf(im_l)
    # 直接计算对数似然和，避免数值不稳定
    log_likelihood = stats.binom.logpmf(no_c, no_a, theoretical_prob)
    return -np.sum(log_likelihood)

# 轻微损伤拟合：设置合理初始值+约束
im_log = np.log(PGA_Values)
neg_ll_slight = partial(neg_log_likelihood_sum, im_l=im_log, no_a=slight_damage_analyses, no_c=slight_damage_collapses)
# L-BFGS-B支持参数下界，sigma>0，beta>0，用1e-6避免0值
res_slight = optimize.minimize(neg_ll_slight, x0=(0.2, 0.5), method='L-BFGS-B', bounds=((1e-6, None), (1e-6, None)))

# 严重损伤拟合
neg_ll_heavily = partial(neg_log_likelihood_sum, im_l=im_log, no_a=heavily_damage_analyses, no_c=heavily_damage_collapses)
res_heavily = optimize.minimize(neg_ll_heavily, x0=(0.4, 0.5), method='L-BFGS-B', bounds=((1e-6, None), (1e-6, None)))

# 倒塌损伤拟合
neg_ll_collapsed = partial(neg_log_likelihood_sum, im_l=im_log, no_a=collapsed_damage_analyses, no_c=collapsed_damage_collapses)
res_collapsed = optimize.minimize(neg_ll_collapsed, x0=(0.5, 0.5), method='L-BFGS-B', bounds=((1e-6, None), (1e-6, None)))

# 打印结果
print("轻微损伤拟合结果：")
print(res_slight)
print("\n严重损伤拟合结果：")
print(res_heavily)
print("\n倒塌损伤拟合结果：")
print(res_collapsed)

# 绘制 fragility 曲线
x = np.linspace(0.1, 0.8, 100)  # 避免log(0)的问题
x_log = np.log(x)

# 轻微损伤曲线
y_slight = stats.norm(loc=np.log(res_slight.x[0]), scale=res_slight.x[1]).cdf(x_log)
# 严重损伤曲线
y_heavily = stats.norm(loc=np.log(res_heavily.x[0]), scale=res_heavily.x[1]).cdf(x_log)
# 倒塌损伤曲线
y_collapsed = stats.norm(loc=np.log(res_collapsed.x[0]), scale=res_collapsed.x[1]).cdf(x_log)

plt.figure(figsize=(10,6))
plt.plot(x, y_slight, label='轻微损伤', color='blue')
plt.plot(x, y_heavily, label='严重损伤', color='orange')
plt.plot(x, y_collapsed, label='倒塌', color='red')

# 绘制观测数据点
plt.scatter(PGA_Values, slight_damage_collapses / slight_damage_analyses, color='blue', alpha=0.5, marker='o')
plt.scatter(PGA_Values, heavily_damage_collapses / heavily_damage_analyses, color='orange', alpha=0.5, marker='s')
plt.scatter(PGA_Values, collapsed_damage_collapses / collapsed_damage_analyses, color='red', alpha=0.5, marker='^')

plt.xlabel('峰值地面加速度(PGA)')
plt.ylabel('损伤概率')
plt.title('损伤易损性曲线')
plt.legend()
plt.grid(True)
plt.show()

关键改进说明 #

• 数值稳定性：通过直接计算对数似然，避免了PMF为0时的数值异常，移除了错误的偏移量操作。
• 参数约束：利用L-BFGS-B的边界约束功能，确保sigma和beta始终为正数，避免不合理的参数值。
• 初始参数优化：根据损伤等级的物理规律设置初始值，帮助优化算法快速收敛到最优解。
• 绘图优化：调整x轴范围避免log(0)错误，使用不同标记区分观测点，提升可视化效果。

内容的提问来源于stack exchange，提问作者Orhun Kalyoncu