Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
sp = ss.gamma.fit(data_shifted, floc=0, optimizer='powell') # scipy's answer is used as an initial guess. Scipy is only correct when there is no censored data
guess = [sp[2], sp[0], self.gamma]
self.initial_guess = guess
k = len(guess)
n = len(all_data)
delta_BIC = 1
BIC_array = [1000000]
runs = 0
bnds = [(0, None), (0, None), (0, min(all_data) - offset)] # bounds on the solution. Helps a lot with stability
while delta_BIC > 0.001 and runs < 5: # exits after BIC convergence or 5 iterations
runs += 1
result = minimize(value_and_grad(Fit_Gamma_3P.LL), guess, args=(failures, right_censored), jac=True, method='L-BFGS-B', bounds=bnds)
params = result.x
guess = [params[0], params[1], params[2]]
LL2 = 2 * Fit_Gamma_3P.LL(guess, failures, right_censored)
BIC_array.append(np.log(n) * k + LL2)
delta_BIC = abs(BIC_array[-1] - BIC_array[-2])
if result.success is True:
params = result.x
self.success = True
self.alpha = params[0]
self.beta = params[1]
self.gamma = params[2]
else:
self.success = False
print('WARNING: Fitting using Autograd FAILED for Gamma_3P. The fit from Scipy was used instead so the results may not be accurate.')
sp = ss.gamma.fit(all_data, optimizer='powell')
self.alpha = sp[2]
self.beta = sp[0]
self.gamma = sp[1]
label = str('Fitted Gamma_2P (α=' + str(round_to_decimals(alpha, dec)) + ', β=' + str(round_to_decimals(beta, dec)) + ')')
if 'color' in kwargs:
color = kwargs.pop('color')
data_color = color
else:
color = 'red'
data_color = 'k'
plt.xlabel('Time')
elif fit_gamma is True:
if __fitted_dist_params is not None:
alpha = __fitted_dist_params.alpha
beta = __fitted_dist_params.beta
gamma = __fitted_dist_params.gamma
else:
from reliability.Fitters import Fit_Gamma_3P
fit = Fit_Gamma_3P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
alpha = fit.alpha
beta = fit.beta
gamma = fit.gamma
gf = Gamma_Distribution(alpha=alpha, beta=beta).CDF(show_plot=False, xvals=xvals)
if 'label' in kwargs:
label = kwargs.pop('label')
else:
label = str('Fitted Gamma_3P\n(α=' + str(round_to_decimals(alpha, dec)) + ', β=' + str(round_to_decimals(beta, dec)) + ', γ=' + str(round_to_decimals(gamma, dec)) + ')')
if 'color' in kwargs:
color = kwargs.pop('color')
data_color = color
else:
color = 'red'
data_color = 'k'
plt.xlabel('Time - gamma')
failures = failures - gamma
def LL(params, T_f, T_rc): # log likelihood function (3 parameter Gamma)
LL_f = 0
LL_rc = 0
LL_f += Fit_Gamma_3P.logf(T_f, params[0], params[1], params[2]).sum() # failure times
LL_rc += Fit_Gamma_3P.logR(T_rc, params[0], params[1], params[2]).sum() # right censored times
return -(LL_f + LL_rc)
# Kaplan-Meier estimate of quantiles. Used in P-P plot.
d = sorted(self._all_data) # sorting the failure data is necessary for plotting quantiles in order
nonparametric = KaplanMeier(failures=failures, right_censored=right_censored, print_results=False, show_plot=False)
self._nonparametric_CDF = 1 - np.array(nonparametric.KM) # change SF into CDF
# parametric models
self.__Weibull_3P_params = Fit_Weibull_3P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
self.Weibull_3P_alpha = self.__Weibull_3P_params.alpha
self.Weibull_3P_beta = self.__Weibull_3P_params.beta
self.Weibull_3P_gamma = self.__Weibull_3P_params.gamma
self.Weibull_3P_BIC = self.__Weibull_3P_params.BIC
self.Weibull_3P_AICc = self.__Weibull_3P_params.AICc
self._parametric_CDF_Weibull_3P = self.__Weibull_3P_params.distribution.CDF(xvals=d, show_plot=False)
self.__Gamma_3P_params = Fit_Gamma_3P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
self.Gamma_3P_alpha = self.__Gamma_3P_params.alpha
self.Gamma_3P_beta = self.__Gamma_3P_params.beta
self.Gamma_3P_gamma = self.__Gamma_3P_params.gamma
self.Gamma_3P_BIC = self.__Gamma_3P_params.BIC
self.Gamma_3P_AICc = self.__Gamma_3P_params.AICc
self._parametric_CDF_Gamma_3P = self.__Gamma_3P_params.distribution.CDF(xvals=d, show_plot=False)
self.__Expon_2P_params = Fit_Expon_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
self.Expon_2P_lambda = self.__Expon_2P_params.Lambda
self.Expon_2P_gamma = self.__Expon_2P_params.gamma
self.Expon_2P_BIC = self.__Expon_2P_params.BIC
self.Expon_2P_AICc = self.__Expon_2P_params.AICc
self._parametric_CDF_Exponential_2P = self.__Expon_2P_params.distribution.CDF(xvals=d, show_plot=False)
self.__Lognormal_3P_params = Fit_Lognormal_3P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
self.Lognormal_3P_mu = self.__Lognormal_3P_params.mu
self.alpha = sp[2]
self.beta = sp[0]
self.gamma = sp[1]
params = [self.alpha, self.beta, self.gamma]
self.loglik2 = LL2
if n - k - 1 > 0:
self.AICc = 2 * k + LL2 + (2 * k ** 2 + 2 * k) / (n - k - 1)
else:
self.AICc = 'Insufficient data'
self.BIC = np.log(n) * k + LL2
self.distribution = Gamma_Distribution(alpha=self.alpha, beta=self.beta, gamma=self.gamma)
# confidence interval estimates of parameters
Z = -ss.norm.ppf((1 - CI) / 2)
hessian_matrix = hessian(Fit_Gamma_3P.LL)(np.array(tuple(params)), np.array(tuple(failures)), np.array(tuple(right_censored)))
covariance_matrix = np.linalg.inv(hessian_matrix)
self.alpha_SE = abs(covariance_matrix[0][0]) ** 0.5
self.beta_SE = abs(covariance_matrix[1][1]) ** 0.5
self.gamma_SE = abs(covariance_matrix[2][2]) ** 0.5
self.alpha_upper = self.alpha * (np.exp(Z * (self.alpha_SE / self.alpha)))
self.alpha_lower = self.alpha * (np.exp(-Z * (self.alpha_SE / self.alpha)))
self.beta_upper = self.beta * (np.exp(Z * (self.beta_SE / self.beta)))
self.beta_lower = self.beta * (np.exp(-Z * (self.beta_SE / self.beta)))
self.gamma_upper = self.gamma * (np.exp(Z * (self.gamma_SE / self.gamma))) # here we assume gamma can only be positive as there are bounds placed on it in the optimizer.
self.gamma_lower = self.gamma * (np.exp(-Z * (self.gamma_SE / self.gamma)))
Data = {'Parameter': ['Alpha', 'Beta', 'Gamma'],
'Point Estimate': [self.alpha, self.beta, self.gamma],
'Standard Error': [self.alpha_SE, self.beta_SE, self.gamma_SE],
'Lower CI': [self.alpha_lower, self.beta_lower, self.gamma_lower],
'Upper CI': [self.alpha_upper, self.beta_upper, self.gamma_upper]}
def LL(params, T_f, T_rc): # log likelihood function (3 parameter Gamma)
LL_f = 0
LL_rc = 0
LL_f += Fit_Gamma_3P.logf(T_f, params[0], params[1], params[2]).sum() # failure times
LL_rc += Fit_Gamma_3P.logR(T_rc, params[0], params[1], params[2]).sum() # right censored times
return -(LL_f + LL_rc)