Reparaciones mínimas.
Modelo NHPP.

El proceso no homogéneo de Poisson

Suponemos un sistema que se mantiene en funcionamiento hasta que se produce un fallo; a continuación se repara de forma que quede en el mismo estado en el que estaba justo antes del fallo, luego se vuelve a poner en servicio hasta el próximo fallo, repitiéndose este procedimiento de forma indefinida. Este modelo se denomina de reparación mínima o ABAO (As Bad As Old).

El modelo formal consiste en una generalización del proceso HPP y se denomina Proceso no homogéneo de Poisson o NHPP (Non Homogeneous Poisson Process), el cual se define en los siguietes términos:

La familia de variables aleatorias \(\{N(t), t\gt 0\}\) es un proceso NHPP si se cumplen las siguientes condiciones:

• \(N(0)=0\)
• \(N(t) \in \mathbb{N}, \forall t\)
• \(N(t_1) \leq N(t_2)\) si \(t_1 \lt t_2\)
• \([N(t_1)-N(0)], [N(t_2)-N(t_1)], [N(t_3)-N(t_2)], \ldots\) son v.a. independientes para \(0 \lt t_1 \lt t_2 \lt t_3, \ldots\), (condición de los incrementos independientes).
• Para cualesquiera \(t \gt s \geq 0\) y \(c \gt 0\) las v.a. \([N(t)-N(s)]\) y \([N(t+c)-N(s+c)]\) tienen la misma distribución de probabilidad, (condición de estacionariedad).
• \(\lim_{h \rightarrow 0} \frac{\Pr(N(t+h)-N(t)=1)}{h}=\lambda(t)\), lo cual se interpreta como que la probabilidad de fallo en un intervalo pequeño es \(\lambda(t) h\).
• \(\lim_{h \rightarrow 0} \frac{\Pr(N(t+h)-N(t) \geq 2)}{h}=0\), lo que se interpreta como que la probabilidad de dos o más fallos en un intervalo pequeño es 0.

El aspecto fundamental que caracteriza al proceso NHPP frente al HPP es que la probabilidad de fallar en un intervalo pequeño es función del tiempo a través de \(\lambda(t)\), la cual que se denomina tasa de fallo o función de intensidad, en contraposición al caso homogéneo, que es constante.

Como en el caso homogéneo, \(N(t)\) es el número de fallos que ha experimentado el sistema durante el intervalo temporal \((0, t]\) y experimentará un incremento de una unidad cada vez que se produzca un fallo. Si el sistema falla en los instantes \(s_1, s_2, \ldots\), estos puntos son las realizaciones de cierto proceso \(\{S_i, i \in \mathbb{N}\}\). Entre dos fallos consecutivos transcurre un tiempo aleatorio definido por \(T_i=T_{i}-T_{i-1}\), lo que define otro proceso estocástico de interés, \(\{T_i, i \in \mathbb{N}\}\).

Se centra el interés en las siguientes variables aleatorias:

• Número de fallos en un intervalo de tiempo \((v_1,v_2]\): se trata de la variable \(N(v_2)-N(v_1)\), cuya distribución es de Poisson de parámetro igual \(\int_{v_1}^{v_2} \lambda(u)du\).
• Tiempo que transcurre entre dos fallos consecutivos: las variables aleatorias \(T_i\) ya no están todas igualmente distribuidas. Si el sistema ha experimentado un fallo en el instante \(v_1\), la función de densidad del tiempo de espera hasta que se produzca el siguiente fallo es \(f_{v_1}(t)=\lambda(v_1+t) e^{-\int_{v_1}^{v_1+t} \lambda(u)du}\), con \(t \gt 0\).

Función de logverosimilitud

Nos planteamos el siguiente esquema muestral. Empezando en \(t=0\), observamos el sistema y vamos anotando los tiempos en los que se producen hasta n fallos, \(s_1, s_2, \ldots, s_n\), instante en el que detenemos las observaciones. Cuando se produce un fallo, detenemos el reloj y procedemos a una reparación mínima, terminada la cual volvemos a poner en marcha el reloj y el sistema.

La función de verosimilitud la definimos como \[ \begin{array}{ccl} L(\lambda(t), s_i) & = & \prod_{i=1}^n \lambda(s_{i-1}+t_i) e^{-\int_{s_{i-1}}^{s_{i-1}+t_i} \lambda(u) du}\\ & = & \prod_{i=1}^n \lambda(s_i) e^{-\int_{s_{i-1}}^{s_i} \lambda(u) du} \\ & = & \prod_{i=1}^n \lambda(s_i) \cdot e^{-\int_0^{s_n} \lambda(u) du}, \end{array} \] donde definimos \(s_0=0\), el momento inicial de la observación. Tomamos logaritmos y obtenemos la función de logverosimilitud, \[ \begin{array}{ccl} l(\lambda(t), s_i) & = & \sum_{i=1}^n \log \lambda(s_i) - \sum_{i=1}^n \int_{s_{i-1}}^{s_i} \lambda(u) du \\ & = & \sum_{i=1}^n \log \lambda(s_i) - \int_0^{s_n} \lambda(u) du. \end{array} \] Recordamos que bajo este patrón de muestreo, \(s_n\) es justo el momento en el que detenemos las observaciones.

El siguiente paso es fijar una expresión para la función de intensidad y estimar los parámetros que contenga.

Modelo NHPP Power-Law

En este caso, \(\lambda(t)=\lambda \beta t^{\beta-1}\), con ambos parámetros \(\lambda, \beta \gt 0\). Calculamos la función de verosimilitud con Maxima.

assume(s[n]>0)$
w(t):=lambda*beta*t^(beta-1)$
L: prod(w(s[i]),i,1,n)*exp(-integrate(w(u),u,0,s[n]));

\[ e^ {- s_{n}^{\beta}\,\lambda }\,\prod_{i=1}^{n}{\beta\,s_{i}^{\beta-1}\,\lambda} \]

En el cálculo anterior, Maxima preguntó si b-1 es positivo o negativo; cualquiera de las dos respuestas llevará al mismo resultado. Ahora es necesario aplicar el logaritmo para obtener la función de logverosimilitud,

logexpand: super$
lv: log(L);

\[ \sum_{i=1}^{n}{\left(\log \lambda+\left(\beta-1\right)\,\log s_{i}+\log \beta\right)}-s_{n}^{\beta}\,\lambda \]

Para obtener el estimador \(\hat{\lambda}\), derivamos respecto de λ, igualamos a cero y despejamos,

mv1: solve(diff(lv,lambda)=0,lambda);

\[ \left[ \lambda={{n}\over{s_{n}^{\beta}}} \right] \]

Hacemos lo propio con β para calcular el estimador \(\hat{\beta}\),

simpsum: true$
mv2: solve(diff(lv,beta)=0,beta);

\[ \left[ \beta={{n}\over{s_{n}^{\beta}\,\log s_{n}\,\lambda-\sum_{i=1}^{n}{\log s_{i}}}} \right] \]

Vemos que β no está completamente aislada porque Maxima no ha sido capaz de gestionar el exponencial; esto se puede resolver sustituyendo λ por el valor antes calculado y así le ayudamos un poco,

mv2: solve(subst(mv1, mv2),beta);

\[ \left[ \beta={{n}\over{s_{n}^{\beta}\,\log s_{n}\,\lambda-\sum_{i=1}^{n}{\log s_{i}}}} \right] \]

Ya tenemos los estimadores de máxima verosimilitud para λ y β.

Ejemplo numérico

Vamos a ajustar el modelo NHPP-Power Law a datos simulados. En primer lugar guardamos la secuencia de los instantes en los que se produjeron los fallos en una variable,

s: [5.66, 8.56, 161.33, 186.29, 1122.23,
    4730.97, 5829.90, 7210.37, 11833.61, 
    14555.25, 19089.61, 20963.65, 21377.01] $

mv(d):=
  block(
    [n:length(d), bet, lam],
    bet: n/(n*log(d[n]) - 
         sum(log(d[i]),i,1,n)),
    lam: n/d[n]^bet,
    [lam, bet])$

mv(s);

\[ \left[ 0.27580068656987 , 0.38645923948952 \right] \]

La lista contiene los estimadores de máxima verosimilitud \(\hat{\lambda}\) y \(\hat{\beta}\), en este orden.

© 2011-2016, TecnoStats.