1- Мои входные данные, которые будут использоваться для построения кривой хрупкости, изображены, как вы можете видеть ниже:
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
