import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az
from scipy.special import gamma as Gamma
from survival_mcmc import weibull_rwm
d = pd.read_csv('shock_absorber.csv'); t = d['km'].to_numpy(float); delta = d['failed'].to_numpy(float); logt = np.log(t)
with pm.Model() as m:
    b0 = pm.Normal('b0', 0, 100); logk = pm.Normal('logk', 0, 5); k = pm.math.exp(logk)
    xb = b0*np.ones(len(t))
    pm.Potential('lik', (delta*(xb+logk+(k-1)*logt) - pt.exp(xb)*t**k).sum())
    pm.Deterministic('k', k); pm.Deterministic('eta', pm.math.exp(-b0/k))
    idata = pm.sample(3000, tune=2000, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
po = idata.posterior; k_p = po['k'].values.ravel(); eta_p = po['eta'].values.ravel()
print('PyMC : shape k = %.2f [%.2f, %.2f]   char. life eta = %.0f [%.0f, %.0f] km   r_hat %.3f' %
      (k_p.mean(), np.percentile(k_p,2.5), np.percentile(k_p,97.5), eta_p.mean(), np.percentile(eta_p,2.5), np.percentile(eta_p,97.5), float(az.rhat(idata)['k'])))
g = weibull_rwm(np.ones((len(d),1)), t, delta, R=60000, burn=10000, seed=1, prior_sd=100.0, logk_sd=5.0, s=0.6)['draws']
gk = np.exp(g[:,1]); geta = np.exp(g[:,0])**(-1/gk)
print('\n%-16s %12s %12s' % ('param','from-scratch','PyMC'))
print('%-16s %12.2f %12.2f' % ('shape k', gk.mean(), k_p.mean()))
print('%-16s %12.0f %12.0f' % ('char. life eta', geta.mean(), eta_p.mean()))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [b0, logk]

Sampling 4 chains for 2_000 tune and 3_000 draw iterations (8_000 + 12_000 draws total) took 9 seconds.

PyMC : shape k = 3.24 [2.15, 4.54]   char. life eta = 25671 [21899, 31541] km   r_hat 1.001

param            from-scratch         PyMC
shape k                  3.32         3.24
char. life eta          25552        25671

Reliability — PyMC cross-check (single-sample Weibull, censored)¶

Results¶