Physics Data Fitting

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Model function
def exponential_decay(t, A, tau, C):
    """y = A * exp(-t/tau) + C"""
    return A * np.exp(-t / tau) + C

# Data with uncertainties
t_data = np.array([0.5, 1, 2, 3, 5, 7, 10, 15, 20])
y_data = np.array([9.8, 8.1, 5.5, 3.9, 2.1, 1.3, 0.7, 0.35, 0.22])
y_err = np.array([0.3, 0.2, 0.2, 0.15, 0.1, 0.1, 0.08, 0.05, 0.04])

# Fit with uncertainties
popt, pcov = curve_fit(exponential_decay, t_data, y_data,
                       sigma=y_err, absolute_sigma=True,
                       p0=[10, 3, 0.1])  # initial guess

# Extract parameters with uncertainties
perr = np.sqrt(np.diag(pcov))
A, tau, C = popt
dA, dtau, dC = perr
print(f"A   = {A:.3f} ± {dA:.3f}")
print(f"τ   = {tau:.3f} ± {dtau:.3f}")
print(f"C   = {C:.3f} ± {dC:.3f}")

# Chi-squared
residuals = (y_data - exponential_decay(t_data, *popt)) / y_err
chi2 = np.sum(residuals**2)
ndof = len(t_data) - len(popt)
chi2_red = chi2 / ndof
print(f"χ²/ndof = {chi2:.2f}/{ndof} = {chi2_red:.2f}")

from lmfit import Model, Parameters

model = Model(exponential_decay)
params = Parameters()
params.add('A', value=10, min=0)         # bounded
params.add('tau', value=3, min=0.01)     # positive definite
params.add('C', value=0, min=-1, max=1)  # bounded

result = model.fit(y_data, params, t=t_data,
                   weights=1/y_err)

# Full report
print(result.fit_report())

# Confidence intervals (profile likelihood)
ci = result.conf_interval()
print("\nConfidence Intervals:")
for name in result.params:
    print(f"  {name}: {ci[name]}")

# Plot with confidence band
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')

# Confidence band via parameter covariance
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)

# Residuals
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-likelihood (Gaussian errors)
            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 comparison
    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

# Example: compare exponential vs power law
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)

    # Compute Jacobian numerically
    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

    # Propagate: Σ_f = J Σ_p J^T
    cov_f = J @ cov_matrix @ J.T

    return f0, np.sqrt(np.diag(cov_f))

# Example: compute half-life from decay constant
def half_life(params):
    tau = params[1]  # decay time
    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}")