from lmfit import Model, Parameters
model = Model(exponential_decay)
params = Parameters()
params.add('A', value=10, min=0)
params.add('tau', value=3, min=0.01)
params.add('C', value=0, min=-1, max=1)
result = model.fit(y_data, params, t=t_data,
weights=1/y_err)
print(result.fit_report())
ci = result.conf_interval()
print("\nConfidence Intervals:")
for name in result.params:
print(f" {name}: {ci[name]}")
fig, axes = plt.subplots(2, 1, figsize=(10, 8), gridspec_kw={'height_ratios': [3, 1]})
t_fine = np.linspace(0, 22, 200)
ax = axes[0]
ax.errorbar(t_data, y_data, yerr=y_err, fmt='ko', capsize=3, label='Data')
ax.plot(t_fine, exponential_decay(t_fine, *popt), 'r-', linewidth=2, label='Best fit')
from scipy.stats import t as t_dist
n_boot = 1000
y_samples = np.zeros((n_boot, len(t_fine)))
rng = np.random.default_rng(42)
for i in range(n_boot):
p_sample = rng.multivariate_normal(popt, pcov)
y_samples[i] = exponential_decay(t_fine, *p_sample)
y_lo = np.percentile(y_samples, 2.5, axis=0)
y_hi = np.percentile(y_samples, 97.5, axis=0)
ax.fill_between(t_fine, y_lo, y_hi, alpha=0.2, color='red', label='95% CI')
ax.set_ylabel('y', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax = axes[1]
ax.errorbar(t_data, residuals, yerr=1, fmt='ko', capsize=3)
ax.axhline(0, color='r', linestyle='--')
ax.set_xlabel('t', fontsize=13)
ax.set_ylabel('Residual (σ)', fontsize=13)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('fit_with_residuals.png', dpi=150, bbox_inches='tight')
def model_comparison(models, data_x, data_y, data_err):
"""Compare models using AIC and BIC."""
results = []
for name, func, p0 in models:
try:
popt, pcov = curve_fit(func, data_x, data_y,
sigma=data_err, absolute_sigma=True, p0=p0)
residuals = (data_y - func(data_x, *popt)) / data_err
chi2 = np.sum(residuals**2)
n = len(data_y)
k = len(popt)
ndof = n - k
log_L = -0.5 * chi2 - 0.5 * n * np.log(2*np.pi) - np.sum(np.log(data_err))
aic = 2*k - 2*log_L
bic = k*np.log(n) - 2*log_L
results.append({
'name': name, 'k': k, 'chi2': chi2,
'chi2_red': chi2/ndof, 'AIC': aic, 'BIC': bic
})
except RuntimeError:
results.append({'name': name, 'k': 0, 'chi2': np.inf,
'chi2_red': np.inf, 'AIC': np.inf, 'BIC': np.inf})
print(f"{'Model':<20} {'k':<4} {'χ²/ndof':<10} {'AIC':<10} {'BIC':<10}")
print("-" * 54)
for r in sorted(results, key=lambda x: x['AIC']):
print(f"{r['name']:<20} {r['k']:<4} {r['chi2_red']:<10.3f} {r['AIC']:<10.2f} {r['BIC']:<10.2f}")
return results
def power_law(t, A, n, C):
return A * t**(-n) + C
models = [
('Exponential', exponential_decay, [10, 3, 0.1]),
('Power law', power_law, [10, 1, 0.1]),
]
model_comparison(models, t_data, y_data, y_err)
def error_propagation(func, params, cov_matrix, *args):
"""
Propagate parameter uncertainties through a function.
Uses the Jacobian: σ²_f = J^T · Σ · J
"""
eps = 1e-8
f0 = func(params, *args)
n_params = len(params)
J = np.zeros((np.size(f0), n_params))
for i in range(n_params):
p_up = params.copy()
p_up[i] += eps
J[:, i] = (func(p_up, *args) - f0) / eps
cov_f = J @ cov_matrix @ J.T
return f0, np.sqrt(np.diag(cov_f))
def half_life(params):
tau = params[1]
return np.array([tau * np.log(2)])
t_half, dt_half = error_propagation(half_life, popt, pcov)
print(f"Half-life: {t_half[0]:.3f} ± {dt_half[0]:.3f}")
