Model Keandalan Perangkat Lunak

# ================================================================
# S11409 - Sesi 3: Fitting Software Reliability Growth Model (SRGM)
# Model: Goel-Okumoto (G-O) NHPP
# Dosen: Riadi Marta Dinata, S.Ti., M.Kom. | ISTN Jakarta
# ================================================================

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import kstest
import warnings
warnings.filterwarnings('ignore')

# ── DATA KEGAGALAN KUMULATIF (bulan, cumulative failures) ──
t_data = np.array([1, 2, 3, 4, 5, 6, 7, 8])
m_data = np.array([5, 9, 12, 14, 15, 15.5, 15.8, 16])  # kumulatif

# ── DEFINISI MODEL ──
def goel_okumoto(t, a, b):
    """Goel-Okumoto Mean Value Function: μ(t) = a(1 - e^(-bt))"""
    return a * (1 - np.exp(-b * t))

def failure_intensity_go(t, a, b):
    """Failure Intensity: λ(t) = ab·e^(-bt)"""
    return a * b * np.exp(-b * t)

def jelinski_moranda_rate(i, phi, N0):
    """J-M Failure Rate: λ_i = φ(N0 - (i-1))"""
    return phi * (N0 - (i - 1))

# ── FITTING PARAMETER MODEL G-O ──
print("="*55)
print("  FITTING GOEL-OKUMOTO MODEL")
print("="*55)

popt, pcov = curve_fit(goel_okumoto, t_data, m_data,
                        p0=[20, 0.5],     # initial guess: a=20, b=0.5
                        bounds=([0, 0], [np.inf, np.inf]),
                        maxfev=10000)
a_fit, b_fit = popt
perr = np.sqrt(np.diag(pcov))

print(f"Parameter a (total bugs): {a_fit:.4f} ± {perr[0]:.4f}")
print(f"Parameter b (detect rate): {b_fit:.4f} ± {perr[1]:.4f}")

# ── PREDIKSI & EVALUASI ──
t_pred = np.linspace(0, 12, 200)  # prediksi hingga bulan 12
m_pred = goel_okumoto(t_pred, a_fit, b_fit)
lam_pred = failure_intensity_go(t_pred, a_fit, b_fit)

# Hitung residuals dan goodness-of-fit
m_fitted = goel_okumoto(t_data, a_fit, b_fit)
residuals = m_data - m_fitted
ss_res = np.sum(residuals**2)
ss_tot = np.sum((m_data - np.mean(m_data))**2)
r_squared = 1 - (ss_res / ss_tot)

print(f"\nGoodness of Fit R²        : {r_squared:.6f}")
print(f"RMSE                      : {np.sqrt(ss_res/len(m_data)):.4f}")

# ── PREDIKSI BUG TERSISA ──
bugs_found = m_data[-1]                       # sampai saat ini
bugs_remaining = a_fit - bugs_found
t_target = 9                                  # target bulan ke-9
m_at_target = goel_okumoto(t_target, a_fit, b_fit)
lam_at_target = failure_intensity_go(t_target, a_fit, b_fit)

print(f"\n--- PREDIKSI ---")
print(f"Total bug diperkirakan      : {a_fit:.1f}")
print(f"Bug ditemukan s/d bulan 8   : {bugs_found:.1f}")
print(f"Bug tersisa (perkiraan)     : {bugs_remaining:.2f}")
print(f"Failure intensity bulan 9   : {lam_at_target:.4f}/bulan")

# ── KAPAN LAYAK RILIS? ──
target_intensity = 0.1  # target failure rate yang bisa diterima
# Cari t ketika λ(t) ≤ target: t = -ln(target / (a*b)) / b
if a_fit * b_fit > target_intensity:
    t_release = -np.log(target_intensity / (a_fit * b_fit)) / b_fit
    print(f"Estimasi waktu siap rilis   : Bulan ke-{t_release:.1f}")
    print(f"  (saat failure intensity ≤ {target_intensity}/bulan)")

# ── VISUALISASI ──
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle('Software Reliability Growth Model - Goel-Okumoto',
             fontsize=14, fontweight='bold')

# Plot 1: Cumulative Failures (Data vs Model)
axes[0].scatter(t_data, m_data, color='#C8102E', s=80,
                zorder=5, label='Data Aktual')
axes[0].plot(t_pred, m_pred, color='#38BDF8', lw=2,
             label=f'G-O Fit (a={a_fit:.1f}, b={b_fit:.2f})')
axes[0].axhline(a_fit, color='#D4A017', linestyle='--',
                label=f'Total bugs ≈ {a_fit:.1f}')
axes[0].set_xlabel('Waktu Pengujian (bulan)')
axes[0].set_ylabel('Kumulatif Kegagalan')
axes[0].set_title('Cumulative Failure vs G-O Model')
axes[0].legend(fontsize=8); axes[0].grid(alpha=0.3)

# Plot 2: Failure Intensity (menurun = reliability meningkat)
axes[1].plot(t_pred, lam_pred, color='#22C55E', lw=2)
axes[1].axhline(target_intensity, color='#F97316',
                linestyle='--', label=f'Target λ = {target_intensity}')
axes[1].fill_between(t_pred, lam_pred, target_intensity,
                      where=(lam_pred <= target_intensity),
                      alpha=0.2, color='#22C55E', label='Zona Aman Rilis')
axes[1].set_xlabel('Waktu Pengujian (bulan)')
axes[1].set_ylabel('Failure Intensity λ(t)')
axes[1].set_title('Failure Intensity (Reliability Growth)')
axes[1].legend(fontsize=8); axes[1].grid(alpha=0.3)

# Plot 3: Residuals
axes[2].bar(t_data, residuals,
           color=['#22C55E' if r >= 0 else '#C8102E' for r in residuals],
           alpha=0.8)
axes[2].axhline(0, color='white', lw=1)
axes[2].set_xlabel('Waktu Pengujian (bulan)')
axes[2].set_ylabel('Residual (Aktual - Prediksi)')
axes[2].set_title(f'Residuals (R² = {r_squared:.4f})')
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('srgm_goel_okumoto.png', dpi=150, bbox_inches='tight')
print("\n✅ Grafik disimpan: srgm_goel_okumoto.png")