Я написал сценарий для подгонки логнормальной функции кумулятивного распределения (CDF) к наблюдаемым данным с использованием MLE, но я столкнулся с проблемами оптимизации. В частности:
• For one damage level, the optimization terminates with unrealistic parameter estimates and doesn’t converge properly.
• For another damage level, the results are logical, and the optimization is successful.
• For the third level, the optimization stops because the maximum number of iterations or function evaluations is exceeded, and the parameter estimates are extremely high or unrealistic.
Результаты показывают, что алгоритм оптимизации испытывает затруднения, возможно, из-за неправильных оценок начальных параметров, числовой нестабильности или проблем с определением функции правдоподобия.
Как правильно подогнать логнормальное распределение к моим данным с помощью MLE, чтобы избежать этих проблем? Существуют ли в Python стандартные методы или встроенные функции, которые могут помочь упростить этот процесс?
1- Мои входные данные, которые будут использоваться для построения кривой хрупкости, изображены, как вы можете видеть ниже: p>
2- Проблема с моим скриптом возникла, когда я его запустил (я попробовал три кривые хрупкости (легкая, сильная и разрушенная)): значения среднего и стандартного отклонения для незначительного и разрушенного нелогичны.
print(res_slight)
final_simplex: (array([[1., 1.],
[1., 1.],
[1., 1.]]), array([inf, inf, inf]))
fun: inf
message: 'Maximum number of function evaluations has been exceeded.'
nfev: 403
nit: 101
status: 1
success: False
x: array([1., 1.])
print(res_heavily)
final_simplex: (array([[0.63531661, 0.88823413],
[0.6353701 , 0.88832485],
[0.63533704, 0.88817963]]), array([203.91066709, 203.91066948, 203.91068102]))
fun: 203.91066708806784
message: 'Optimization terminated successfully.'
nfev: 60
nit: 32
status: 0
success: True
x: array([0.63531661, 0.88823413])
print(res_collapsed)
final_simplex: (array([[8.67945396e+20, 3.53694860e+01],
[7.84728875e+20, 3.53008691e+01],
[7.54836798e+20, 3.52584876e+01]]), array([183.64705108, 183.64785074, 183.6479341 ]))
fun: 183.64705108393747
message: 'Maximum number of iterations has been exceeded.'
nfev: 1823
nit: 1000
status: 2
success: False
x: array([8.67945396e+20, 3.53694860e+01])
Вот почему мои кривые хрупкости выглядят так, как показано на рисунке:
[img]https://i. sstatic.net/fNiNdX6t.png[/img]
мой код для получения кривых хрупкости:
from functools import partial
import numpy as np
from scipy import optimize, stats
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])
im_log = np.log(PGA_Values)
#Slight Damage Data Results
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 Data Results
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 Data Results
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
theoretical_fragility_function = stats.norm(np.log(sigma), beta).cdf(im_l)
likelihood = stats.binom.pmf(no_c, no_a, theoretical_fragility_function)
log_likelihood = np.log(likelihood)
log_likelihood_sum = np.sum(log_likelihood + 1e+10)
print(FORMAT_STRING.format(sigma, beta, log_likelihood_sum))
return -log_likelihood_sum
# Slight Damage Data fitting using MLE
neg_log_likelihood_sum_partial_slight = partial(neg_log_likelihood_sum, im_l=np.log(PGA_Values), no_a=slight_damage_analyses, no_c=slight_damage_collapses)
res_slight = optimize.minimize(neg_log_likelihood_sum_partial_slight, (1, 1), method="Nelder-Mead")
# Heavily Damage Data fitting using MLE
neg_log_likelihood_sum_partial_heavily = partial(neg_log_likelihood_sum, im_l=np.log(PGA_Values), no_a=heavily_damage_analyses, no_c=heavily_damage_collapses)
res_heavily = optimize.minimize(neg_log_likelihood_sum_partial_heavily, (1, 1), method="Nelder-Mead")
# Collapsed Data fitting using MLE
neg_log_likelihood_sum_partial_collapsed = partial(neg_log_likelihood_sum, im_l=np.log(PGA_Values), no_a=collapsed_damage_analyses, no_c=collapsed_damage_collapses)
res_collapsed = optimize.minimize(neg_log_likelihood_sum_partial_collapsed, (1, 1), method="Nelder-Mead")
# Plotting Fragility Curves
import matplotlib.pyplot as plt
x = np.linspace(0, 1.0, 100)
# Slight Damage
y_slight = stats.norm(np.log(res_slight.x[0]), res_slight.x[1]).cdf(np.log(x))
plt.plot(x, y_slight, label='Slight Damage')
# Heavily Damage
y_heavily = stats.norm(np.log(res_heavily.x[0]), res_heavily.x[1]).cdf(np.log(x))
plt.plot(x, y_heavily, label='Heavily Damage')
# Collapsed
y_collapsed = stats.norm(np.log(res_collapsed.x[0]), res_collapsed.x[1]).cdf(np.log(x))
plt.plot(x, y_collapsed, label='Collapsed')
plt.scatter(PGA_Values, slight_damage_collapses / slight_damage_analyses, color='blue', alpha=0.5)
plt.scatter(PGA_Values, heavily_damage_collapses / heavily_damage_analyses, color='orange', alpha=0.5)
plt.scatter(PGA_Values, collapsed_damage_collapses / collapsed_damage_analyses, color='r', alpha=0.5)
Подробнее здесь: https://stackoverflow.com/questions/790 ... ood-estima
