Stochastic processes

Renewal process

import matplotlib.pyplot as plt
import numpy as np
from relife.lifetime_model import Weibull
from relife.stochastic_process import RenewalProcess

distrib = Weibull(5, 0.03)
renewal_process = RenewalProcess(distrib)

rf_timeline, rf = renewal_process.renewal_function(150, nb_steps=5000)
print(rf_timeline.shape, rf.shape)

rd_timeline, rd = renewal_process.renewal_density(150, nb_steps=5000)
print(rd_timeline.shape, rd.shape)

(5000,) (5000,)
(5000,) (5000,)

fig, ax = plt.subplots(ncols=2, nrows=1,  figsize=(10, 5))

ax[0].set_title("Renewal function")
ax[0].set_xlabel("t")
ax[0].plot(rf_timeline, rf)

ax[1].set_title("Renewal density")
ax[1].set_xlabel("t")
ax[1].plot(rd_timeline, rd)

plt.show()

tf = 150
nb_samples = 100

renewal_sample = renewal_process.sample(tf=150, size=100, seed=10)
print(renewal_sample.select(sample_id=0).sample_id)
print(renewal_sample.select(sample_id=0).timeline)
print(renewal_sample.select(sample_id=0).time)
print(renewal_sample.select(sample_id=0).event)

[0 0 0 0 0]
[ 25.45401812  62.5772016  100.35485328 125.64522842 150.        ]
[25.45401812 37.12318348 37.77765168 25.29037514 24.35477158]
[ True  True  True  True False]

timeline, mean_nb_events = renewal_sample.mean_nb_events()
plt.plot(rf_timeline, rf, label="renewal function")
plt.plot(timeline, mean_nb_events, alpha=0.6, label="mean nb of events (sample size = 100)")
plt.legend()

<matplotlib.legend.Legend at 0x780d7009d7d0>

Computation of equilibrium age distribution. Over a sufficiently long period of time, the varariance of life spans “spreads” failures over time, resulting in a stabilization of the rate of occurrence of failures over time. When this stationary regime is reached, the working population has an age distribution that no longer varies over time: this is the equilibrium age distribution.

from relife.lifetime_model import EquilibriumDistribution

eq_distrib = EquilibriumDistribution(distrib)

fig, ax = plt.subplots()
ax.set_ylim(top=0.08)
distrib.plot.pdf(ax=ax, end_time=50, label="Lifetimes distribution")
eq_distrib.plot.pdf(ax=ax, end_time=50, label="Equilibrium lifetimes distribution")
plt.show()

renewal_sample = renewal_process.sample(tf=150, size=5000, seed=21)
sample_age_eq = [renewal_sample.select(sample_id=i).time[-1] for i in range(1000)]
sample_age = [renewal_sample.select(sample_id=i).time[-2] for i in range(1000)]

t = np.linspace(0, 50, num=1000)

fig, ax = plt.subplots()
ax.set_ylim(top=0.08)

ax.plot(t, eq_distrib.pdf(t), label="Real equilibrium lifetimes distribution")
ax.hist(sample_age_eq, density=True, alpha=0.5, color=ax.lines[-1].get_color(), label="Sample distribution of equilibrium lifetimes")

ax.plot(t, distrib.pdf(t), label="Real equilibrium lifetimes distribution")
ax.hist(sample_age, density=True, alpha=0.5, color=ax.lines[-1].get_color(), label="Sample distribution of equilibrium lifetimes")

ax.set_xlabel("Lifetimes")
ax.set_ylabel("Density")
ax.legend(loc="upper left")
plt.show()