Análisis de la Influencia de los Datos Censurados en el Cálculo del Intervalo Óptimo entre Mantenimientos Preventivos y Estudio Comparativo de Varios Métodos de Estimación

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Máster Máster Universitario en Ingeniería Aeronáutica Análisis de la Influencia de los Datos Censurados en el Cálculo del Intervalo Óptimo entre Mantenimientos Preventivos y Estudio Comparativo de Varios Métodos de Estimación Autor: Antonio Á. Cuenca Ávalos Tutor: Antonio J. Sánchez Herguedas Dpto. Organización Industrial y Gestión de Empresas I Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Máster Máster Universitario en Ingeniería Aeronáutica Análisis de la Influencia de los Datos Censurados en el Cálculo del Intervalo Óptimo entre Mantenimientos Preventivos y Estudio Comparativo de Varios Métodos de Estimación Autor: Antonio Á. Cuenca Ávalos Tutor: Antonio J. Sánchez Herguedas Profesor Contratado Doctor Dpto. Organización Industrial y Gestión de Empresas I Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Máster: Análisis de la Influencia de los Datos Censurados en el Cálculo del Intervalo Óptimo entre Mantenimientos Preventivos y Estudio Comparativo de Varios Métodos de Estimación Autor: Antonio Á. Cuenca Ávalos Tutor: Antonio J. Sánchez Herguedas El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos En primer lugar, me gustaría agradecerle a mi tutor, Antonio J. Sánchez Herguedas, el haberme acompañado en este último tramo de mis estudios de máster, prestándome su ayuda cuando la necesitaba y descubriéndome todas estas herramientas matemáticas tan útiles que pueden materializarse en aplicaciones reales. Una vez más, querría dale las gracias a mis padres, a mi hermano y al resto de mi familia, por todo el apoyo mostrado durante toda mi vida, en especial, en los momentos más duros. Gracias a ellos, he llegado a donde estoy y tengo un futuro prometedor por delante. Finalmente, no puedo cerrar este apartado sin darle las gracias a todos mis amigos. Tanto a los de Córdoba, como a los de Cádiz, quienes siguen demostrando lo que es la amistad verdadera, a pesar de habernos ido a vivir a diferentes ciudades, y hacen parecer que no ha pasado el tiempo cuando nos volvemos a reunir. Además, no podían faltar los agradecimientos a mis nuevos amigos de Sevilla, que no solo han demostrado un gran compañerismo, ayudándome en todo momento, si no que le han dado el broche de oro a mi época universitaria. A todos vosotros os dedico mi trabajo. Antonio Ángel Cuenca Ávalos Sevilla, 2022 I Resumen Una de las decisiones más importantes que debe tomar un ingeniero es cuándo realizar los mantenimientos preventivos de un equipo, con el fin de que no fallen debido al desgaste. Históricamente, se han definido los intervalos entre mantenimientos preventivos a partir de los datos de fallo recogidos de la observación de varios equipos en funcionamiento. No obstante, una estimación más precisa puede realizarse teniendo en cuenta los datos censurados producidos en esa misma observación. El objetivo de este proyecto es estudiar cómo influye la inclusión de dichos datos censurados en el cálculo del intervalo óptimo entre mantenimientos preventivos. Para ello, se hace uso de un método extraído de la literatura existente, en el que se supone al equipo en un modelo de tres estados: operativo, correctivo y preventivo. En el primero se obtienen beneficios, en los otros dos, se producen costes. La fórmula que calcula el valor del intervalo óptimo entre mantenimientos preventivos hace uso tanto de estos retornos, como de los parámetros de forma, escala y umbral resultantes de ajustar la muestra de datos de fallo a una distribución Weibull. Para obtener los parámetros de Weibull, se emplean tres métodos de estimación diferentes: la estimación de máxima verosimilitud, el método de Newton-Raphson y la regresión del rango de medianas. El primero de ellos es uno de los métodos con mejor desempeño según la bibliografía revisada. El método de Newton- Raphson, por su parte, es una versión del método anterior que permite la convergencia de la función de verosimilitud a una buena velocidad en las situaciones en las que la estimación de máxima verosimilitud no funciona correctamente. Por último, la regresión del rango de medianas es uno de los métodos más utilizados históricamente por los ingenieros, debido a su sencillez. A lo largo de este proyecto, se incorporan los resultados obtenidos por cada uno de estos métodos y se establecen comparaciones. Por último, cabe destacar que se incluye un capítulo en el que se recogen varios métodos que han sido mencionados a lo largo del documento y se explican brevemente, haciendo hincapié en los tres métodos de estimación usados para calcular los parámetros de Weibull, con el fin de facilitar la comprensión del lector. III Índice de Tablas 1.1 Ventajas del algoritmo EM 4 1.2 Desventajas del algoritmo EM 5 4.1 Parámetros comunes a todos los casos de estudio 22 4.2 Datos de fallo para el caso de desgaste 22 4.3 Datos censurados para el caso de desgaste 23 4.4 Parámetros de Weibull obtenidos mediante el MLE para el caso de desgaste con censura 23 4.5 Resultados del estudio de velocidad de convergencia para el desgaste con datos censurados 25 4.6 Resumen de los resultados obtenidos para el caso de desgaste con datos censurados 25 4.7 Parámetros de Weibull obtenidos mediante el MLE para el caso de desgaste sin censura 27 4.8 Incremento, número de orden ajustado y RMide los datos de fallo sin tener en cuenta la censura 28 4.9 Incremento, número de orden ajustado y RMide los datos de fallo teniendo en cuenta la censura 29 4.10 Resultados del estudio de velocidad de convergencia para el desgaste sin datos censurados 32 4.11 Resumen de los resultados obtenidos para el caso de desgaste sin datos censurados 32 4.12 Valores finales de τ034 4.13 Datos de fallo para el caso de fallos aleatorios 34 4.14 Datos censurados para el caso de fallos aleatorios 34 4.15 Resumen de los resultados obtenidos para el caso de fallos aleatorios 36 4.16 Datos de fallo para el caso de fallos infantiles 37 4.17 Datos censurados para el caso de fallos infantiles 37 4.18 Resumen de los resultados obtenidos para el caso de fallos infantiles 38 XI Notación A Tiempo medio que el sistema permanece en el estado operativo antes de que se produzca el fallo B Tiempo medio que el sistema permanece en el estado correctivo C Tiempo medio que el sistema permanece en el estado preventivo Cn iCoeficiente binomial empleado en la estimación L-moment dPosición anterior a la posición ren la estimación L-moment F(τ τ τ)Matriz de tiempo de permanencia F(t)Función de distribución acumulativa de fallos del equipo f(t)Función de densidad de probabilidad de fallos del equipo G(tc) Función de distribución acumulativa del tiempo en el que el equipo se encuentra bajo mantenimiento correctivo g(tc) Función de densidad de probabilidad del tiempo en el que el equipo se encuentra bajo mantenimiento correctivo h(θ) Derivada de la función de verosimilitud respecto a un estimador cualquiera h′(θ) Segunda derivada de la función de verosimilitud respecto a un estimador cualquiera J(tp) Función de distribución acumulativa del tiempo en el que el equipo se encuentra bajo el mantenimiento preventivo j(tp) Función de densidad de probabilidad del tiempo en el que el equipo se encuentra bajo el mantenimiento preventivo IMatriz unidad L(α,β)Función de verosimilitud lnL(α,β)Función de verosimilitud logarítmica mNúmero de transiciones m′ rMuestra correspondiente a un momento TL m0 1Muestra simplificada para una posición 1 m0 2Muestra simplificada para una posición 2 nTotalidad de observaciones tomadas, tamaño de la muestra N1,N2 N3,N4Parámetros de cálculo de h(θ)yh′(θ) N5 PMatriz de probabilidad de transición entre estados pij Probabilidad ligada a un valor total de los retornos producidos en un estado i y durante su transición al estado j Qi(1)Retorno en un solo paso Qi(m)Retorno acumulado desde el estado ien mtransiciones qiNúmero de orden ajustado de un elemento i qi-1 Número de orden del dato de fallo anterior al elemento i XIII XIV Notación R1Ingresos obtenidos durante el estado operativo R2Costes derivados del estado correctivo R3Costes derivados del estado preventivo Rij Gastos producidos de la transición de un estado i a un estado j RMi Función de distribución observada, rango de la mediana del elemento i rPosición determinada en la estimación L-moment rij Valor total de los retornos producidos en un estado i y durante su transición al estado j tiValor de la observación i V(m)Forma matricial del retorno medio acumulado vi(m)Retorno medio acumulado Xi:n Estadístico iésimo dentro del tamaño nde la muestra αFactor de forma de Weibull ˆ αLM Estimador del factor de forma según la estimación L-moment βFactor de escala de Weibull ˆ βLM Estimador del factor de escala según la estimación L-moment Γ(·)Función Gamma γFactor de localización o umbral de Weibull δiIndicador de censura εError impuesto θkValor inicial de un estimador genérico θk+1 Valor de un estimador genérico en la siguiente iteración µ′ rMomento TL π(θ)Distribución a priori τIntervalo entre mantenimientos preventivos τ0Intervalo óptimo entre mantenimientos preventivos 1 Estudio de la literatura existente En este capítulo, se realiza una revisión de la literatura existente acerca del trato de los datos censurados en sus diferentes tipos, a la hora de realizar estudios de fiabilidad. El objetivo es analizar los diferentes métodos que han utilizado varios expertos a la hora de incluir en sus análisis la información que ofrecen los datos censurados. 1.1 Tipos de datos censurados Los datos censurados hacen referencia a aquella información de la que no se dispone debido a que se desconoce cuando ocurre cierto evento clave. Particularizando en los análisis de fiabilidad, puede realizarse un estudio sobre el tiempo de fallo de una muestra de todos los activos de una planta. De dicho estudio, puede extraerse información de relevancia sobre cuándo pueden fallar dichos activos, lo cual es importante, por ejemplo, a la hora de establecer los intervalos de tiempo entre mantenimientos preventivos. La censura ocurre cuando no se puede conocer el momento exacto del fallo individual de todos los activos de la planta [ 47 ]. Existen diferentes tipos de censura [12]: • Censura por la izquierda. También conocida como left censoring, en inglés. Este tipo de censura indica que el evento de relevancia ha ocurrido en algún instante antes de que el activo sea objeto de estudio, desconociéndose el momento exacto [20]. • Censura por la derecha. Por el contrario, cuando el evento relevante ocurre en un momento desconocido, posterior a cierto limite de tiempo definido, se dice que existe una censura por la derecha, o right censoring, como se conoce en inglés [ 47 ]. Por su parte, existen varios subtipos de censura por la derecha: – Censura tipo 1. En este caso, existe un tiempo fijo de censura, que puede ser el tiempo de estudio [ 19 ]. De este modo, el evento de relevancia quedará determinado en todos los sujetos de estudio que lo sufran durante el análisis. Todos aquellos que no lo experimenten en el tiempo fijado, puede que lo hagan tiempo después de que termine el estudio, aunque será indeterminado, quedando dicho dato censurado. – Censura tipo 2. Por su parte, la censura tipo 2 aparece en aquellos análisis que definen su duración según la proporción de activos que sufren el fallo [ 19 ]. Por ejemplo, si en cierto análisis se estudian seis equipos, se puede establecer que el estudio termina cuando cuatro de ellos registren el fallo, quedándose indeterminados los tiempos de fallo correspondientes a los otros dos equipos. Es por ello que en la literatura puede encontrarse textos que se refieren a este tipo de censura como “censura de fallo” (failure censoring) [ 47 ]. La figura 1.1 ilustra el ejemplo propuesto. – Censura tipo 3. Mientras que los dos tipos anteriores de censura se dan en aquellos análisis en los que los activos que son objeto de estudio son introducidos al principio del periodo experimental, este tercer tipo de censura tiene lugar en aquellos en los que cada activo presenta un instante inicial diferente. No obstante, también presentan un tiempo fijo de censura. Es por ello que este tipo de censura recibe el nombre de “censura tipo 1 generalizada” [19]. 1 2Capítulo 1. Estudio de la literatura existente • Censura de intervalo. Cuando en un mismo estudio se dan los dos tipos de censuras anteriores, se dice que existe una censura en intervalo (interval censoring) [20]. Por otro lado, hay que tener en cuenta que no son lo mismo los datos censurados que los datos truncados. En el truncamiento, solo se tienen en cuenta aquellos activos cuyo evento ocurre dentro de una ventana de observación. De este modo, no se dispone de información de todo activo cuyo evento quede fuera de dicha ventana. Esto es lo que contrasta con la censura, en la que hay al menos información parcial sobre cada activo de estudio [ 20 ]. Este documento se centra en los datos censurados, pero se han descrito los datos truncados ya que pueden ser mencionados en adelante. 0 1 2 3 4 5 6 0 5 10 15 20 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 35 40 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F (a) Censura tipo 1. 0 1 2 3 4 5 6 0 5 10 15 20 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 35 40 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F (b) Censura tipo 2. 0 1 2 3 4 5 6 0 5 10 15 20 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F 0 1 2 3 4 5 6 0 5 10 15 20 25 30 35 40 Tiempo de estudio (semanas) Fin del estudio Casos de estudio A B C D E F (c) Censura tipo 3. Figura 1.1 Representación gráfica de los tipos de censura por la derecha. 1.2 Revisión de la literatura 3 1.2 Revisión de la literatura Una vez definidos los diferentes tipos de datos censurados, se hace uso de la revisión de la literatura realizada para exponer las diferentes técnicas que han empleado numerosos autores a la hora de incluir los datos censurados en sus estudios, como se indica al principio del capítulo. La distribución Weibull es un modelo de probabilidad ampliamente utilizado tanto en los estudios de vida útil, como en la teoría de fiabilidad [ 61 ]. Gran parte de la literatura revisada hace uso de ella para tratar las muestras con datos censurados. Cabe recordar que la distribución de Weibull depende de tres parámetros: el de forma ( α ), el de escala o vida característica ( β ) y el de localización ( γ ), también conocido como valor umbral o vida garantizada. Cuando este último tiene un valor nulo, se dispone de una distribución Weibull de dos parámetros [25]. El trabajo realizado por Teimouri et al. [ 26 ] establece una comparativa entre diferentes métodos de estimación propuestos para resolver una distribución Weibull. Estos son la estimación de máxima verosimilitud, el método del momento logarítmico, el método del percentil, el método de los momentos y el método de estimación L-moment, propuesto por los autores. El estudio numérico consiste en un conjunto de 10 000 individuos con los que se hacen varios estudios con muestras de diferente tamaño. El modelo se ajusta a una distribución Weibull de dos parámetros y se utiliza el sesgo y el error de mínimos cuadrados como criterios de comparación. Finalmente, los autores concluyen que la estimación de máxima verosimilitud es la que mejor comportamiento presenta, para muestras lo suficientemente grandes. Le sigue la estimación L-moment, que actúa correctamente en la mayoría de casos estudiados. Se remarca el hecho de que, en ciertos casos, estos dos métodos alcanzan los mismos valores del error de mínimos cuadrados. El resto de métodos no alcanzan tan buenos resultados. Según Genschen et al. [ 60 ], históricamente, el método de los mínimos cuadrados ordinarios ha sido una alternativa empleada por los ingenieros a la estimación de máxima verosimilitud, debido a su fácil programación. El método más común a la hora de implementar el enfoque de mínimos cuadrados ordinarios es la regresión del rango de medianas. En su artículo, establecen una comparativa entre ambos métodos para poner de manifiesto las ventajas que tiene cada uno de ellos. Antes de empezar el estudio, los autores parten de las siguientes ideas: • Según la experiencia de los autores, es difícil encontrar un método de estimación que supere el buen funcionamiento del de máxima verosimilitud. • Con el paso del tiempo, los métodos de mínimos cuadrados ordinarios han ido dejando de estar incluidos en los softwares comerciales más comunes, a diferencia del método de máxima verosimilitud. Por ello, su programación a día de hoy no es mucho más sencilla que la de la estimación de máxima verosimilitud. • Tanto los estimadores de máxima verosimilitud, como los de la regresión del rango de medianas están sesgados. En el segundo caso, el sesgo suele estar dominado por la varianza; por lo que, cuando se compara la precisión global de los estimadores de ambos métodos, suelen dar mejores resultados los de máxima verosimilitud. • La regresión del rango de medianas dibuja la línea atravesando los puntos trazados. El método de máxima verosimilitud lo hace cuando corresponde, pero la regresión del rango de medianas lo hace incluso cuando no debe. Eso puede llevar a resultados erróneos. • La estimación de máxima verosimilitud presenta propiedades óptimas en muestras lo suficientemente grandes. • La estimación de máxima verosimilitud puede aplicarse en situaciones de censura más complicada, como la de intervalo o truncamiento. Como criterios de comparación, los autores emplean los parámetros de forma y de escala, así como varios cuantiles de Weibull. Se realizan estudios con muestras de datos con censura tipo 1 y 2. Finalmente, los autores llegan a la conclusión de que la estimación de máxima verosimilitud presenta un mejor comportamiento que la estimación por regresión del rango de medianas a lo largo de la mayoría del estudio. Referente a la censura tipo 4Capítulo 1. Estudio de la literatura existente 1, la estimación de máxima verosimilitud tiene ciertas ventajas, ya que no ignora la información censurada. Igualmente, posee una mayor precisión que la regresión del rango de medianas ante datos censurados del tipo 2, aunque las diferencias son menores que las detectadas en las simulaciones con datos censurados del tipo 1. Por su parte, Shen et al. proponen en su artículo [ 5 ] un nuevo método para llevar a cabo los estudios de fiabilidad con datos censurados. Dicho método hace uso de una serie de funciones de regresión, analizada mediante un enfoque Bayesiano y un algoritmo de Monte Carlo basado en Cadenas de Markov (MCMC, por sus siglas en inglés), para transformar los datos censurados en datos de fallo. Los autores aseguran que, de este modo, se simplifica el estudio de fiabilidad, manteniendo la precisión computacional. Este nuevo método tendría aplicación en los casos donde la muestra es reducida y la censura es alta. Por otro lado, el uso del algoritmo de esperanza-maximización para resolver la estimación de máxima verosimilitud puede traer problemas de convergencia si el tamaño de la muestra no es lo suficientemente grande, lo cual es consistente con lo indicado en [ 60 ]. Igualmente, el algoritmo de Monte Carlo basado en Cadenas de Markov puede suponer una alternativa al algoritmo de esperanza-maximización si los cálculos se vuelven demasiado complejos y no se tiene una precisión adecuada. Otros autores que comparan varios métodos de estimación de parámetros para una distribución Weibull son Avijit Joarder et al. En su artículo [ 6 ], indican que los estimadores de máxima verosimilitud no pueden obtenerse de forma explícita, haciendo uso de un proceso iterativo, por lo que proponen la llamada estimación de máxima verosimilitud aproximada, la cual presenta forma explícita en sus estimadores y, además, puede ser facilmente computada. La contribución principal de este artículo es establecer una comparativa entre estas dos estimaciones y la estimación bayesiana. La comparación se lleva a cabo mediante numerosas simulaciones de Monte Carlo y se utilizan el sesgo y el error de mínimos cuadrados como criterio comparativo, al igual que en [26]. Como conclusión general, llegan a que todos los métodos funcionan correctamente menos cuando la muestra y el tiempo son demasiados pequeños. Por su parte, la estimación de máxima verosimilitud aproximada se comporta de forma muy similar a la estimación de máxima verosimilitud original, aunque presenta intervalos de confianza ligeramente más cortos. Por otro lado, la estimación bayesiana registra peores comportamientos, además, para muestras muy grandes. Finalmente, indican que si no se tienen ninguna información de los parámetros desconocidos, la estimación de máxima verosimilitud aproximada puede usarse para la mayoría de casos prácticos. Los artículos expuestos hasta el momento coinciden en que el método que mejor comportemiento presenta, de forma generalizada, es el de máxima verosimilitud. Siguiendo este mismo hilo, el artículo escrito por Ferreira y Silva [ 56 ] tiene como objetivo validar la aplicación del algoritmo de esperanza-maximización como método para resolver la estimación de máxima verosimilitud de un conjunto de datos incompletos, censurados por la derecha. Además de lograr dicho objetivo, los autores incluyen una lista de ventajas y desventajas del algoritmo de esperanza-maximización que puede resultar interesante a la hora de decidir el algoritmo a utilizar en la parte práctica de este proyecto. Dicha lista se encuentra plasmada en las tablas 1.1 y 1.2 Tabla 1.1 Ventajas del algoritmo EM. Ventajas Converge en un gran número de condiciones, excepto si no se hace una buena elección de los valores iniciales del parámetro a maximizar. El trabajo analítico requerido es más sencillo que el de otros métodos. Es relativamente sencillo de programar e implementar. Permite controlar la convergencia y los errores de programación durante las iteracciones. 1.2 Revisión de la literatura 5 Tabla 1.2 Desventajas del algoritmo EM. Desventajas Puede converger lentamente si la información sufre una fuerte censura o truncamiento. No garantiza la convergencia global si existen varios puntos máximos. La estimación obtenida depende de la solución inicial. Por otro lado, en el artículo desarrollado por Balakrishnan y Mitra [ 51 ], se dispone de ciertos datos truncados por la izquierda y censurados por la derecha. Los autores hacen uso de un algoritmo EM Gradiente para ajustar la distribución Weibull a dicho conjunto de datos. La razón por la que usan dicho algoritmo está relacionada con la primera desventaja plasmada en la tabla 1.2, ya que suele emplearse para acelerar la convergencia del modelo [ 13 ]. Además, con fines comparativos, se realiza otro ajuste mediante el método de Newton-Raphson, el cual también acelera la convergencia de la estimación de máxima verosimilitud y es más sencillo. En dicho artículo, los autores llegaron a la conclusión de que el algoritmo EM Gradiente presenta algunos problemas de convergencia para muestras de alrededor de 100 individuos e inferiores, lo cual coincide con lo indicado por Shen et al. [ 5 ]. Igualmente, se detectan problemas del mismo tipo si se toman valores iniciales de los parámetros muy alejados de los reales. Por último, queda comentar que, aunque se alcanzan resultados muy similares con ambos métodos, se registran pequeñas diferencias en los intervalos de confianza y en las probabilidades de confianza. Mientras que los intervalos de confianza para el parámetro de localización, son similares para los dos algoritmos, en el caso del parámetro de escala, el algoritmo EM Gradiente presenta intervalos de confianza más amplios que el método de Newton-Raphson. Del mismo modo, para mayores porcentajes de truncamiento, mayores probabilidades de cobertura del parámetro de localización se registran. Finalmente, cabe comentar que Robert Ross [ 38 ] [ 39 ] y Hideo Hirose [ 17 ], proponen en sus trabajos un método para reducir el sesgo de los estimadores de máxima verosimilitud, para muestras pequeñas, de forma que mejora el comportamiento de la estimación. Además, Zhang et al. [ 22 ] proponen una modificación de los métodos propuestos por Ross y Hirose aplicada a la estimación de mínimos cuadrados. Respecto a los datos con censura del tipo 2, Riffat Jabeen et al. [ 36 ] realizan un estudio en el que comparan la eficiencia de dos métodos con los que obtener los estimadores de la muestra: el BLUE (Best Linear Unbiased Estimator) y el método Gupta. Finalmente, llegan a la conclusión de que el método de Gupta puede reemplazarse por el BLUE, ya que se alcanzan valores de eficiencia relativa bastante alta. Esto es útil cuando se conocen las esperanzas, pero la matriz de varianza-covarianza es desconocida. Por otro lado, como James K. Starling et al. comentan en su artículo [ 19 ], en muchas aplicaciones reales suele darse una censura de tipo 1 generalizada. Son conscientes de los problemas de sesgo que presenta la estimación de máxima verosimilitud y proponen un método mejorado combinando dicha estimación con un método de sobremuestreo. La conclusión del estudio es que los métodos propuestos reducen el sesgo en los estimadores del parámetro de forma, pero aumentan el sesgo en el de escala. No obstante, mediante una reducción en la divergencia de Kullback–Leibler, se puede solucionar el problema de aumento de sesgo del estimador del parámetro de escala. En ciertas situaciones, los datos de la muestra pueden mostrar cierta heterogeneidad. Esto significa que la fiabilidad de los individuos de la muestra depende de parámetros muy diversos. Para tratar dicha heterogeneidad, se suele usar la distribución de mezcla Weibull de dos componentes. De este modo, para afrontar los problemas derivados de una alta censura en los datos, Ducros y Pamphile [ 57 ] proponen un método de bootstrap bayesiano llamado Maximización de la Restauración Bayesiana (BRM, por sus siglas en inglés). Las ventajas de este método son que presenta mejores resultados en términos de precisión respecto al algoritmo EM y a su variante, el algoritmo EM estocástico. Además, al no ser un método iterativo, requiere de tiempos menores de computación. Por otro lado, Musleh y Helu, en su artículo [ 37 ], indican que la distribución Weibull no tiene un buen desempeño en aquellas muestras que disponen de tasas de fallo no monótonas, tal y como puede ser la unimodal. En estos casos, la distribución Weibull inversa es un modelo más adecuado. En dicho artículo, 12 Capítulo 2. Descripción de los métodos . Figura 2.1 Ejemplo de gráfica obtenida con los datos facilitados en [33]. Los parámetros α y β se pueden obtener fácilmente de la gráfica. Mientras que el primero comparte valor con la pendiente de la línea de tendencia, el factor de escala se calcula mediante la exponencial de la división entre la ordenada en el origen (b) y la pendiente de la línea de tendencia [4]. β=exp−b α(2.17) Por su parte, para calcular el valor de γ es necesario incluir otra línea de tendencia, en este caso, cuadrática. Para obtener el valor adecuado de gamma, hay que incluirlo en la expresión del eje X, de forma que ahora se vea representado por ln(ti−γ) , indicando ti las horas de funcionamiento de cierto elemento i antes de producirse el fallo. El γ óptimo es aquel que consiga hacer 0 el término cuadrático de la expresión de la nueva línea de tendencia cuadrática. Cabe indicar que para poder alcanzar dicho valor de γ , es necesario que el término cuadrático de la expresión sea negativo [4]. 2.6 Método de Monte Carlo basado en cadenas de Markov Existen multitud de métodos de Monte Carlo. En particular el método de Monte Carlo basado en cadenas de Markov es uno de los más usados en la estadística bayesiana, con el fin de obtener muestras de distribuciones de probabilidad específicas [35]. Este método presenta un buen desempeño, respecto al método de Monte Carlo ordinario, ante funciones más complicadas gracias al uso de las cadenas de Markov. Para poder generar dichas cadenas, se debe hacer uso de un algortimo [35]. El más popular es el de Metrópolis-Hastings [14]. Una vez que las cadenas converjan, se dispondrá de una nueva muestra, ajustada a la función de densidad de la distribución uniforme correspondiente [14], y se podrán realizar los cálculos deseados [35]. 2.7 Método BLUE Este método parte de la regresión de mínimos cuadrados ordinarios [ 44 ]. Se habla del método BLUE cuando los estimadores obtenidos por mínimos cuadrados cumplen el teorema de Gauss-Markov, que indica que cuando un estimador presente determinadas condiciones (homocedasticidad, combinación lineal de variables observadas y no observadas, esperanza condicionada de las perturbaciones nula...) es el estimador lineal e 2.8 Método Gupta 13 insesgado óptimo [ 50 ] y presenta una varianza mínima [ 27 ]. Precisamente, BLUE son las siglas en inglés para Best Linear Unbiased Estimators. 2.8 Método Gupta Este método fue propuesto por A. K. Gupta en su trabajo [ 16 ], donde presentó unos estimadores lineales alternativos. El método BLUE requiere que los valores exactos de la media, varianza y covarianza sean conocidos. Esto dificulta su uso en aquellas situaciones en las que se complique el cálculo de la covarianza o en las que muestra sea demasiado grande, dificultando la inversión de la matriz de varianza-covarianza [ 48 ]. En estas situaciones se puede recurrir al método de Gupta. Los estimadores lineales alternativos se obtienen mediante la sustitución de la matriz de varianza-covarianza por una matriz unidad, I [ 36 ], facilitando el cálculo de los estimadores cuando no se disponga de los parámetros requeridos para obtener los estimadores lineales e insesgados óptimos. 2.9 Estimación de Bayes Existen una gran variedad de métodos bayesianos, empleados por los analistas en los estudios de fiablilidad [29]. Esta estimación consiste en combinar los datos reales de una muestra con las incertidumbres propias de las distribuciones a priori que se tengan en cuenta, con el fin de obtener nuevas distribuciones de parámetros conocidas como a posteriori, de las que se obtiene los estimadores. Esto provoca que la validez de los resultados esté ligada tanto a la validez del modelo, como de la distribución a priori [10] [20] [29]. 2.10 Maximización de la restauración bayesiana La maximización de la restauración bayesiana es un método desarrollado por M. Bacha y G. Celeux [ 52 ] basado en el concepto del bootstrap propuesto por B. Efron [11]. Este método es capaz de sobrellevar fuertes censuras en la muestra, ya que se calculan los estimadores a partir de una nueva muestra obtenida de la simulación de la muestra original y la distribución a priori [ 57 ]. En eso radica la ventaja de este método, en que tiene en cuenta la información a priori desde el primer paso [52]. Su funcionamiento es sencillo, el algoritmo realiza tres pasos [52]: •Toma una muestra θde la distribución a priori π(θ). •Obtiene la información que falta, conociendo la muestra observada y el parámetro θobtenido. • Finalmente, obtiene el estimador correspondiente maximizando la verosimilitud de la muestra completa. 3 Antecedentes El objetivo de este proyecto es estudiar cómo influyen los datos censurados en el cálculo del tiempo óptimo entre mantenimientos preventivos, así como en el comportamiento de la distribución Weibull. Esta sección pretende tanto presentar el marco en el que se desarrolla el modelo, como exponer el trasfondo matemático que posee. 3.1 Marco teórico En el supuesto sobre el que se desarrolla el proyecto, los equipos pueden transitar entre tres estados diferentes, que representan el comportamiento del sistema ante el fallo. Estos son operativo, correctivo y preventivo. El estado operativo supone que el equipo está operando en sus condiciones de diseño, produciendo ciertos ingresos. Los otros dos estados se dan dependiendo del momento en el que se dé el fallo del activo. De este modo, si transcurre un intervalo de tiempo predeterminado τ , conocido como intervalo preventivo, sin que se produzca el fallo, el equipo pasa al estado preventivo, caracterizado por los costes de la mano de obra y de las piezas de recambio derivados del mantenimiento preventivo, además de los originados por la inactividad del equipo. Por otro lado, si el fallo tiene lugar antes de alcanzar el tiempo τ , el equipo pasará al estado correctivo, donde quedará inoperativo hasta que sea reparado. Una vez finalizadas las acciones de mantenimiento, el equipo vuelve al estado operacional. Cada vez que el equipo transita entre dos estados, se acumulan los costes o los ingresos en forma de retornos. El cálculo del intervalo preventivo óptimo se lleva a cabo maximizando el retorno medio acumulado. La figura 3.1 representa la transición entre los tres estados expuestos. Estado 1: Operativo Estado 2: Correctivo Estado 3: Preventivo T = τ T < τ Figura 3.1 Transición entre los tres estados del modelo.. Como se ha comentado anteriormente, cada estado supone unos ingresos o unos costes. En este proyecto, se le denomina R1 a los ingresos obtenidos durante el estado operativo. Por su parte, R2 y R3 representan los costes derivados de los estados correctivo y preventivo, respectivamente. Finalmente, cabe comentar que la transición entre estados supone también unos costes. Se representan como Rij los gastos producidos de la transición de un estado ia un estado j. 15 16 Capítulo 3. Antecedentes 3.2 Cálculo del intervalo preventivo Una vez definido el modelo en el que se basa este proyecto, llega el momento de exponer la base matemática que existe detrás del cálculo del parámetro τ. Como se ha comentado anteriormente, un equipo puede oscilar entre tres estados diferentes. En un estado i , y en la posterior transición a un estado j , se producen ciertos costes o ingresos. En un total de m transiciones, se genera un retorno acumulado desde el estado i denominado Qi(m) . Cabe destacar que un sistema evoluciona de forma impredecible desde un estado inicial definido. Por ello, el parámetro de interés a la hora de realizar el desarrollo matemático es el retorno medio acumulado vi(m) , el cual se encuentra representado por la figura 3.2. Dicho parámetro atiende a la siguiente expresión: vi(m) = vi(1)+ n ∑ j=1 vj(m−1)pij (3.1) Se puede observar que la expresión 3.1 dispone de dos términos. El primero representa la transición desde el estado inicial i hasta el estado siguiente, j . El retorno en un solo paso Qi(1) es una variable aleatoria, que puede alcanzar diversos valores ( ri1(1) , ri2(1) , ..., rin(1) ). De este modo, el retorno medio en un solo paso es producto de la multiplicación de dichos valores por su respectivas probabilidades (pi1,pi2, ..., pin). vi(1) = n ∑ j=i rij(1)·pij (3.2) Por su parte, el segundo término de la expresión 3.1 hace referencia al valor del retorno medio acumulado de las m−1 transiciones restantes donde, igualmente, se multiplican los valores aleatorios por sus probabilidades. t Estado 1: Operativo Estado 2: Correctivo Estado 3: Preventivo T < τT = τ Estado 1: Operativo Retorno medio acumulado Transiciones 1 2 3 4 Figura 3.2 Evolución del retorno medio acumulado.. Finalmente, para poder calcular el retorno medio acumulado en m pasos desde cualquier estado i inicial, se puede expresar la ecuación 3.1 de forma matricial: V(m) = V(1)+PV(m−1)(3.3) El parámetro P es la matriz de probabilidad de transición entre estados, la cual atiende a la siguiente expresión: 3.2 Cálculo del intervalo preventivo 17 P=  0F(τ)1−F(τ) 1 0 0 1 0 0  (3.4) Al definir P , aparece otra nueva matriz que definir: la de tiempo de permanencia F . Cada componente Fij representa el tiempo medio que el sistema permanece en el estado i , antes de pasar al j . De este modo, se tiene que: F=  0Aτ B0 0 C0 0 (3.5) De igual modo, en la expresión 3.5 aparecen otros tres términos en los que cabría profundizar. Dichos parámetros atienden a las siguientes expresiones: A=1 F(τ)Zτ 0 t·f(t)dt (3.6) B=Z∞ 0 tc·g(tc)dtc(3.7) C=Z∞ 0 tp·j(tp)dtp(3.8) Como se puede observar, las tres expresiones anteriores dependen de tres funciones de densidad de probabilidad. El parámetro A depende de la función f(t) , que está relacionada con la función de distribución acumulativa F(t) de fallos del equipo. Por otro lado, el parámetro B depende de la función g(tc) que, a su vez, depende del tiempo necesario para finalizar la reparación, el cual está relacionado con la función de distribución G(tc) . Por último, j(tp) es la función de densidad de probabilidad que hace acto de presencia en la expresión 3.8, la cual está relacionada con la función de distribución del tiempo en el que el equipo se encuentra bajo mantenimiento preventivo. Disponiendo de todos estos parámetros, también se puede definir la matriz de retornos R. R=  0r12 r13 r21 0 0 r31 0 0  =  0R1A+R12 R1τ+R13 R2B+R12 0 0 R3C+R12 0 0  (3.9) A continuación, se procede a resolver el sistema de ecuaciones diferencial que supone la expresión 3.3. Para ello, se aplica la transformada Z a ambos lados de la expresión. Por definición, la transformada Z de una secuencia cualquiera x(y), es una función de variable compleja que se define mediante la siguiente serie de Laurent: Z[x(m)] = ∞ ∑ n=0 x(n)z−n;z∈C{0}(3.10) Después de aplicar la transformada Z a la expresión 3.3 y de hacer los cambios pertinentes, se alcanza la siguiente expresión: Z[V(m)] = 1 z−1(I−z−1P)−1V(1)(3.11) Llegados a este punto, para obtener el retorno medio acumulado de m transiciones solo quedaría revertir la transformación. Teniendo en cuenta que el modelo en el que se basa este proyecto oscila entre tres estados, el vector V(1)se puede expresar de la siguiente manera: V(1) =   v1(1) v2(1) v3(1) =  r12F(τ)+r13 [1−F(τ)] v2(1) v3(1) (3.12) Como se observa la evolución del sistema desde el estado operacional, se realizan los cálculos con el primer componente del vector V(1). De este modo, para alcanzar definitivamente el valor óptimo de τ, hay 18 Capítulo 3. Antecedentes que obtener la expresión explícita de v1(1) , derivarla respecto de τ e igualarla a cero. Por un lado, la expresión explícita de τes la siguiente: v1(m) = 1 4(2m+1+(−1)m−1)R1Zτ 0 t f (t)dt+R12F(τ)+(R1τ+R13)[1−F(τ)]+ (2m−1−(−1)m−1)R2Z∞ 0 tcg(tc)dtc+R21F(τ)+R3Z∞ 0 tpj(tp)dtp+R31[1−F(τ)] (3.13) Una vez definida la expresión 3.13, se procede a derivarla respecto de τ: dv1(m) dτ=1 4{f(τ)[(2m+1+(−1)m−1)(R12 −R13)+(2m−1−(−1)m−1) (R2B+R21 −R3C−R31)]}+R1 4[2m+1+(−1)m−1][1−F(τ)] (3.14) Como se ha comentado al principio del capítulo, se hace uso de una distribución Weibull de tres parámetros para estudiar el modelo propuesto. La función de densidad de probabilidad f(t) y la función de distribución acumulativa F(t)atienden a las siguientes expresiones: f(t) = α β"t−γ βα−1#exp−t−γ βα(3.15) F(t) = 1−exp−t−γ βα(3.16) Finalmente, teniendo en cuenta que los parámetros α , β y γ son mayores que 0, igualando la expresión 3.14 a cero y sustituyendo en ella las expresiones 3.15 y 3.16 se alcanza la expresión con la que se obtiene el valor óptimo del intervalo preventivo. (τ0−γ)α−1=βα α·−R1 R12 −R13 +2m−1−(−1)m−1 2m−1+(−1)m−1·(R2B+R21 −R3C−R31) (3.17) Cabe destacar que la expresión 3.17 es independiente del número de transiciones para un valor de m par. No obstante, cuando m presenta un valor impar, el valor de τ0 va descendiendo hasta alcanzar el valor independiente de m cuando m→∞ . De igual modo, para que la expresión tenga una solución válida es necesario que el denominador de la fracción cuyo numerador es −R1 sea negativo. De este modo, α y β han de ser positivos. Este desarrollo matemático ha sido extraído de [1], [2], [3] y [4]. 3.3 Teoría de fallo Como se acaba de indicar, la expresión 3.17 permite calcular el intervalo óptimo entre mantenimientos preventivos. No obstante, el comportamiento de dicha expresión puede verse alterado por la naturaleza de la tasa de fallos que muestren los equipos que conformen el conjunto de datos de estudio [1]. La tasa de fallo se define como la probabilidad de que cierto elemento o equipo falle [ 32 ]. Una de las formas más comunes de representar el comportamiento de dicha tasa a lo largo del tiempo es la conocida curva de bañera, la cual se puede dividir en tres tramos: fallos infantiles, vida útil y desgaste [31]. Como se puede observar en la figura 3.3, la curva de bañera es, en realidad, resultado de la superposición de tres distribuciones de fallos, las cuales se corresponden con las tres regiones indicadas anteriormente [ 43 ]. 3.3 Teoría de fallo 19 Figura 3.3 Curva de bañera [43]. Mientras que una tasa de fallos decreciente denota un aumento de la fiabilidad de los equipos, debido a la sustitución de los elementos más débiles que fallaron demasiado pronto, una tasa de fallos creciente supone que empiezan a predominar los efectos del desgaste en los equipos [ 49 ]. Por su parte, la tasa de fallo constante está relacionada con el desempeño de vida útil del equipo, donde los fallos ocurren de forma aleatoria [ 31 ]. La región de la curva de bañera en la que se encuentre el conjunto de datos influye directamente en los valores de los parámetros de Weibull, de ahí que afecte a los resultados de la expresión 3.17. En particular, es el parámetro de forma el que indica la región de la curva de bañera en la que se encuentra el conjunto de datos. Esto tiene sentido, ya que α está relacionada con la pendiente de la curva. [ 49 ]. De esta forma, se tiene que: •α<1 representa una tasa de fallo decreciente, relacionada con fallos infantiles, de forma que una gran cantidad de fallos se dan al principio del marco temporal. Suelen estar relacionados con defectos de fabricación [ 49 ]. Como se ha comentado anteriormente, la tasa de fallo disminuye al sustituir dichos elementos defectuosos, demostrando los equipos un correcto desempeño tras la sustitución [32]. •α=1 está relacionado con una tasa de fallo constante. De este modo, los fallos que se dan lugar son aleatorios [ 49 ], a menudo relacionados con fluctuaciones fortuitas de tensión que superan la resistencia del componente [43]. •α>1 indica que la tasa de fallo es creciente, de modo que el desgaste en el elemento se está haciendo notar, dando lugar al final de la vida útil del mismo [ 49 ]. De este modo, valores de α superiores a la unidad indican el momento de realizar un mantenimiento preventivo, con el fin de alargar la vida útil del activo. En lo referente al cálculo del intervalo óptimo entre mantenimientos preventivos, cabe destacar que se descarta el caso en el que α=1 , ya que representa el caso en el que los fallos son aleatorios, por lo que no se podría hacer una estimación adecuada. Además, esto se encuentra matemáticamente reflejado en la expresión 3.17; ya que, si α=1 , cualquier número puede ser solución o, por el contrario, ninguno, dependiendo del segundo término de la expresión [1]. Por otro lado, un valor de α inferior a la unidad supone que el máximo del retorno medio acumulado para m transiciones se obtiene para un τ0 de valor infinito. De este modo, se concluye que no procede realizar mantenimientos preventivos en esta región de la curva de bañera [ 1 ]. Esto, además, tiene sentido, ya que, los fallos infantiles están relacionados con defectos de fabricación que provocan el fallo del equipo a las pocas horas de funcionamiento. En estas situaciones, corresponde sustituir el elemento defectuoso, lo cual está relacionado con el mantenimiento correctivo. De este modo, el último caso que queda es aquel en el que α>1 , donde la expresión 3.17 presenta un buen desempeño. De este modo, si se tiene un valor τ inferior al valor óptimo calculado, la función v1(m) 20 Capítulo 3. Antecedentes presenta una tendencia creciente, que alcanza su máximo cuando τ=τ0 y empieza a decrecer cuando τ>τ0 [1]. Este comportamiento se encuentra ilustrado en la figura 3.4, mediante un ejemplo de diez transiciones. Figura 3.4 Comportamiento del retorno medio acumulado para diez transiciones [1]. Por lo tanto, los valores de α superiores a la unidad son los que permiten observar el comportamiento teórico de τ0según el valor de m, el cual fue expuesto al final de la sección 3.2. Figura 3.5 Comportamiento teórico de τ0para diez transiciones [1]. La figura 3.5 sigue con el ejemplo de las diez transiciones presentado en la figura 3.4. Como se puede observar, el valor óptimo de τ es constante para los valores pares de m , mientras que los impares presentan un intervalo óptimo superior, cuyo valor se va aproximando al del número de transiciones pares en cuanto más transiciones se llevan a cabo. Por último, cabe comentar que tiene sentido que este sea el caso en el que se puede calcular el intervalo óptimo para el mantenimiento preventivo, ya que los equipos se encuentran en la región de desgaste y el objetivo del mantenimiento preventivo es poner a punto los equipos antes de que se produzca el fallo por desgaste. A continuación, en el capítulo 4, se incluye una aplicación matemática para cada uno de los tres casos, con el fin de obtener una serie de conclusiones. 4 Aplicación matemática En este capítulo, se desarrollan tres casos de estudio que demuestran la aplicabilidad real del desarrollo matemático y la expresión 3.17, basado tanto en la teoría de fallo, como en el marco teórico expuestos a lo largo del capítulo 3. 4.1 Base común a los diferentes casos de estudio Todos los casos de estudios están ligados al marco expuesto en la sección 3.1. De este modo, se dispone de un conjunto de datos diferente para cada uno de los estudios realizados. Dichos datos agrupan información sobre la operatividad de ciertos equipos que se someten a un determinado periodo de observación, en el que se impone un intervalo entre mantenimientos preventivos de 2 000 h ( τ=2 000 ). Durante la toma de datos, todos los activos empezaron a ser observados en el estado de operación. Si registraban un fallo antes de llegar a las 2 000 h, se enviaban al mantenimiento correctivo. Por el contrario, si alcanzaban el valor de τ impuesto, se enviaban a mantenimiento preventivo, obteniendo así datos censurados por la derecha, ya que se desconoce cuándo podría haberse producido el fallo. Otro caso en el que pueden obtenerse datos censurados es aquel en el que el equipo estuviese funcionando hasta el fin del tiempo de análisis, sin registrar ningún fallo ni alcanzar las 2 000 h de funcionamiento. Los tres conjuntos de datos presentan 83 datos de fallo y 38 censurados. El objetivo es obtener el intervalo óptimo entre mantenimientos preventivos y, para ello, se ajustan los conjuntos de datos a una distribución Weibull, tanto de dos, como de tres parámetros. Como se indicaba en la expresión 3.17, el valor de τ0 depende de los parámetros de Weibull. Por ello, se hace uso de tres métodos matemáticos para encontrar el valor de α , β y, en el caso de tres parámetros, γ . Estos son la Estimación de Máxima Verosimilitud, el algoritmo de Newton-Raphson y el Método del Rango de Medianas. La selección de métodos realizada se debe a lo expuesto a lo largo de la sección 1.2. De la revisión de la literatura se extrajo que la estimación de máxima verosimilitud es uno de los métodos más extendidos y de los que mejor comportamineto se registra en general. Por su parte, el método de Newton-Raphson presenta un método iterativo que mejora el comportamiento de la estimación de máxima verosimilitud cuando la función de verosimilitud logarítmica no puede maximizarse de forma analítica. Por último, el Método de Rango de Medianas ha sido históricamente muy utilizado por los ingenieros debido a su fácil programación. Todo esto es por lo que se considera interesante emplear estos tres métodos para obtener los parámetros de Weibull, comparando desempeños y resultados entre ellos. Ademas de las variables de Weibull, para obtener los intervalos óptimos es necesario definir aquellos parámetros relacionados con el retorno de cada transición y estado, así como los parámetros B y C relacionados con los mantenimientos (correctivo y preventivo) de los equipos. Todos ellos fueron presentados a lo largo de la sección 3.2 y sus valores se encuentran reflejados en la tabla 4.1. Finalmente, se analizan los resultados obtenidos. 21 28 Capítulo 4. Aplicación matemática Tabla 4.8 Incremento, número de orden ajustado y RMide los datos de fallo sin tener en cuenta la censura. Orden Horas de funcionamiento Tipo Incremento Nºde orden RMi 1 190 F 1 1 0.011400651 2 276 F 1 2 0.027687296 3 296 F 1 3 0.043973941 4 409 F 1 4 0.060260586 5 429 F 1 5 0.076547231 6 430 F 1 6 0.092833876 7 437 F 1 7 0.109120521 8 454 F 1 8 0.125407166 9 481 F 1 9 0.141693811 10 492 F 1 10 0.157980456 11 498 F 1 11 0.174267101 12 499 F 1 12 0.190553746 13 543 F 1 13 0.206840391 14 552 F 1 14 0.223127036 15 552 F 1 15 0.239413681 16 552 F 1 16 0.255700326 17 552 F 1 17 0.271986971 18 577 F 1 18 0.288273616 19 603 F 1 19 0.304560261 20 604 F 1 20 0.320846906 21 604 F 1 21 0.33713355 22 612 F 1 22 0.353420195 23 619 F 1 23 0.36970684 24 658 F 1 24 0.385993485 25 675 F 1 25 0.40228013 26 683 F 1 26 0.418566775 27 696 F 1 27 0.43485342 28 702 F 1 28 0.451140065 29 742 F 1 29 0.46742671 30 754 F 1 30 0.483713355 31 773 F 1 31 0.5 32 797 F 1 32 0.516286645 33 812 F 1 33 0.53257329 34 836 F 1 34 0.548859935 35 881 F 1 35 0.56514658 36 889 F 1 36 0.581433225 37 912 F 1 37 0.59771987 38 913 F 1 38 0.614006515 39 942 F 1 39 0.63029316 40 974 F 1 40 0.646579805 41 994 F 1 41 0.66286645 42 994 F 1 42 0.679153094 43 1014 F 1 43 0.695439739 44 1015 F 1 44 0.711726384 45 1024 F 1 45 0.728013029 46 1025 F 1 46 0.744299674 47 1041 F 1 47 0.760586319 48 1105 F 1 48 0.776872964 49 1183 F 1 49 0.793159609 50 1203 F 1 50 0.809446254 51 1211 F 1 51 0.825732899 Continúa en la siguiente página 4.2 Caso de estudio 1: Desgaste 29 Tabla 4.8 – Continúa desde la página anterior Orden Horas de funcionamiento Tipo Incremento Nºde orden RMi 52 1236 F 1 52 0.842019544 53 1238 F 1 53 0.858306189 54 1240 F 1 54 0.874592834 55 1249 F 1 55 0.890879479 56 1274 F 1 56 0.907166124 57 1295 F 1 57 0.923452769 58 1304 F 1 58 0.939739414 59 1312 F 1 59 0.956026059 60 1343 F 1 60 0.972312704 61 1345 F 1 61 0.988599349 Tabla 4.9 Incremento, número de orden ajustado y RMide los datos de fallo teniendo en cuenta la censura. Orden Horas de funcionamiento Tipo Incremento Nºde orden RMi 1 45 R - - - 2 84 R - - - 3 84 R - - - 4 103 R - - - 5 121 R - - - 6 176 R - - - 7 190 F 1.0517 1.0517 0.00619 8 199 R - - - 9 203 R - - - 10 207 R - - - 11 211 R - - - 12 215 R - - - 13 219 R - - - 14 223 R - - - 15 276 F 1.1198 2.1726 0.01541 16 296 F 1.1198 3.2915 0.02464 17 409 F 1.119891443 4.411398467 0.033866544 18 429 F 1.119891443 5.531289911 0.04309135 19 430 F 1.119891443 6.651181354 0.052316156 20 437 F 1.119891443 7.771072797 0.061540962 21 454 F 1.119891443 8.89096424 0.070765768 22 481 F 1.119891443 10.01085568 0.079990574 23 492 F 1.119891443 11.13074713 0.08921538 24 498 F 1.119891443 12.25063857 0.098440186 25 499 F 1.119891443 13.37053001 0.107664992 26 508 R - - - 27 514 R - - - 28 519 R - - - 29 543 F 1.155632659 14.52616267 0.117184207 30 552 F 1.155632659 15.68179533 0.126703421 31 552 F 1.155632659 16.83742799 0.136222636 32 552 F 1.155632659 17.99306065 0.14574185 33 552 F 1.155632659 19.14869331 0.155261065 34 577 F 1.155632659 20.30432597 0.16478028 35 603 F 1.155632659 21.45995863 0.174299494 36 604 F 1.155632659 22.61559129 0.183818709 37 604 F 1.155632659 23.77122395 0.193337924 Continúa en la siguiente página 30 Capítulo 4. Aplicación matemática Tabla 4.9 – Continúa desde la página anterior Orden Horas de funcionamiento Tipo Incremento Nºde orden RMi 38 612 F 1.155632659 24.92685661 0.202857138 39 619 F 1.155632659 26.08248927 0.212376353 40 638 R - - - 41 655 R - - - 42 658 F 1.184166799 27.26665607 0.22213061 43 661 R - - - 44 675 F 1.199156252 28.46581232 0.232008339 45 683 F 1.199156252 29.66496857 0.241886067 46 696 F 1.199156252 30.86412482 0.251763796 47 702 F 1.199156252 32.06328108 0.261641525 48 735 R - - - 49 742 F 1.215361067 33.27864214 0.271652736 50 754 F 1.215361067 34.49400321 0.281663947 51 760 R - - - 52 764 R - - - 53 766 R - - - 54 769 R - - - 55 73 F 1.286852894 35.7808561 0.292264054 56 790 R - - - 57 797 F 1.306350665 37.08720677 0.303024767 58 812 F 1.306350665 38.39355743 0.313785481 59 836 F 1.306350665 39.6999081 0.324546195 60 867 R - - - 61 881 F 1.327420837 41.02732893 0.335480469 62 889 F 1.327420837 42.35474977 0.346414743 63 909 R - - - 64 912 F 1.349919495 43.70466927 0.357534343 65 913 F 1.349919495 45.05458876 0.368653944 66 942 F 1.349919495 46.40450826 0.379773544 67 974 F 1.349919495 47.75442775 0.390893145 68 986 R - - - 69 994 F 1.374918005 49.12934576 0.402218664 70 994 F 1.374918005 50.50426376 0.413544183 71 1011 R - - - 72 1014 F 1.401877181 51.90614094 0.425091771 73 1015 F 1.401877181 53.30801812 0.436639359 74 1024 F 1.401877181 54.70989531 0.448186947 75 1025 F 1.401877181 56.11177249 0.459734534 76 1041 F 1.401877181 57.51364967 0.471282122 77 1105 F 1.401877181 58.91552685 0.48282971 78 1111 R - - - 79 1167 R - - - 80 1183 F 1.467080771 60.38260762 0.494914396 81 1188 R - - - 82 1188 R - - - 83 1203 F 1.540434809 61.92304243 0.507603315 84 1211 F 1.540434809 63.46347724 0.520292234 85 1236 F 1.540434809 65.00391205 0.532981154 86 1238 F 1.540434809 66.54434686 0.545670073 87 1240 F 1.540434809 68.08478167 0.558358992 88 1249 F 1.540434809 69.62521648 0.571047912 89 1274 F 1.540434809 71.16565129 0.583736831 Continúa en la siguiente página 4.2 Caso de estudio 1: Desgaste 31 Tabla 4.9 – Continúa desde la página anterior Orden Horas de funcionamiento Tipo Incremento Nºde orden RMi 90 1295 F 1.540434809 72.7060861 0.59642575 91 1304 F 1.540434809 74.24652091 0.60911467 92 1312 F 1.540434809 75.78695572 0.621803589 93 1343 F 1.540434809 77.32739052 0.634492508 94 1345 F 1.540434809 78.86782533 0.647181428 95 1395 R - - - 96 1396 R - - - 97 1407 F 1.658929795 80.52675513 0.660846418 98 1413 F 1.658929795 82.18568492 0.674511408 99 1421 F 1.658929795 83.84461472 0.688176398 100 1442 F 1.658929795 85.50354451 0.701841388 101 1447 F 1.658929795 87.16247431 0.715506378 102 1458 F 1.658929795 88.8214041 0.729171368 103 1479 F 1.658929795 90.4803339 0.742836358 104 1499 F 1.658929795 92.13926369 0.756501348 105 1505 R - - - 106 1523 F 1.7565139 93.89577759 0.770970161 107 1556 F 1.7565139 95.65229149 0.785438974 108 1577 F 1.7565139 97.40880539 0.799907787 109 1600 F 1.7565139 99.16531929 0.8143766 110 1621 F 1.7565139 100.9218332 0.828845413 111 1656 F 1.7565139 102.6783471 0.843314226 112 1701 F 1.7565139 104.434861 0.85778304 113 1745 F 1.7565139 106.1913749 0.872251853 114 1765 F 1.7565139 107.9478888 0.886720666 115 1781 F 1.7565139 109.7044027 0.901189479 116 1812 F 1.7565139 111.4609166 0.915658292 117 1881 F 1.7565139 113.2174305 0.930127105 118 1923 F 1.7565139 114.9739444 0.944595918 119 1970 F 1.7565139 116.7304583 0.959064731 120 2000 R - - - 121 2000 R - - - A continuación, al igual que con el estudio anterior, se incluye la figura 4.5, que plasma tanto los valores de Weibull obtenidos, así como las gráficas de donde provienen. De nuevo, se puede observar cómo el valor de γhace cero el término cuadrático y que eso disminuye el valor de αyβ. Figura 4.5 Resultados obtenidos mediante MRR para el caso de desgaste sin censura. 32 Capítulo 4. Aplicación matemática Finalmente, se añaden los resultados obtenidos con el método de Newton-Raphson. De nuevo, se realiza el estudio de velocidad de convergencia valiéndose de los valores de γcalculados con los otros dos métodos. Figura 4.6 Iteraciones del método de Newton con γ=113.21, para el caso de desgaste con censura. Figura 4.7 Iteraciones del método de Newton con γ=62.90, para el caso de desgaste con censura. Tabla 4.10 Resultados del estudio de velocidad de convergencia para el desgaste sin datos censurados. Estudio de la velocidad de convergencia γcalculada con MLE MRR Iteración de convergencia 6 6 α2.13 2.27 β1 032.36 1 090.39 γ113.21 62.90 Como se puede observar en las figuras 4.6 y 4.7, el método converge en la sexta iteración en ambos estudios, mucho más rápido que en el caso con datos censurados. No obstante, la diferencia de los otros dos parámetros se hace notar mucho más que en la tabla 4.5. La conclusión es la misma que la de la sección anterior, se elige el valor de γ calculado mediante la estimación de máxima verosimilitud. Además, es remarcable el hecho de que el parámetro N4 tiene un valor nulo, lo cual es congruente con la expresión 4.5, ya que solo se tienen en cuenta los datos censurados en su cálculo. A continuación, se incluye una tabla similar a la tabla 4.6, en la que se plasman los resultados obtenidos por los tres métodos. Igualmente, la figura 4.8 ofrece una comparativa entre las gráficas y los valores alcanzados del intervalo entre mantenimientos preventivos. Tabla 4.11 Resumen de los resultados obtenidos para el caso de desgaste sin datos censurados. Método α β γ MLE 2.13 1 032.36 113.21 Newton-Raphson 2.13 1 032.36 113.21 MRR 2.47 848.12 62.895 4.2 Caso de estudio 1: Desgaste 33 (a) Intervalo óptimo entre mantenimientos preventivos según la estimación de máxima verosimilitud. (b) Intervalo óptimo entre mantenimientos preventivos según el método de Newton-Raphson. (c) Intervalo óptimo entre mantenimientos preventivos según la regresión del rango de medianas. Figura 4.8 Intervalos óptimos entre mantenimientos preventivos para el caso de desgaste con censura. Llegados a este punto, se pueden alcanzar ciertas conclusiones comparando los valores de Weibull plamados en las tablas 4.6 y 4.11. Como se puede observar en ellas, los valores de los tres parámetros son inferiores en el caso en el que no se tiene en cuenta la censura. En particular, destaca la disparidad entre los valores de β. Esto provoca que los valores de τ0sean menores en este caso. Como se puede comprobar en la figura 4.8, el valor óptimo del intervalo entre mantenimientos preventivos ha pasado de estar entre unas 1 009 h y unas 977 h a estar entre unas 840 h y unas 814 h. Si bien escoger unos intervalos de tiempo inferiores forman parte de un enfoque más conservador, ya que cuanto antes se lleven a cabo los mantenimientos menos probable es se produzca el fallo por desgaste, lo óptimo es estirar lo máximo posible el intervalo entre mantenimientos, con el fin de aumentar los beneficios percibidos. Por ello, se prefieren los resultados obtenidos al incluir los datos censurados, ya que se consideran más precisos y los intervalos son mayores. 34 Capítulo 4. Aplicación matemática Por otro lado, en el caso de la regresión de rango de medianas, se puede advertir que la diferencia entre los valores de los parámetros de Weibull es aún mayor, y en cada parámetro es diferente. Mientras que en el estudio con censura el valor de α es de 2.09, en este caso alcanza un valor de 2.47. Por su parte, el valor de β se ve disminuido en casi 340 unidades al despreciar los datos censurados. El parámetro γ también se ve disminuido. Esta enorme discrepacia entre los dos casos se debe a la valiosa información que se pierde al eludir los datos censurados. Aunque el método obtenga sus resultados de los datos de fallo, los datos censurados modifican su rango. Al no disponer de ellos, se pierde una gran cantidad de precisión. De este modo, el hecho de que se le sume mayor imprecisión a un método que de por sí es menos exacto que los otros dos, da lugar a que no se tomen en cuenta sus resultados. Para concluir este estudio, puede decirse que los valores de τ0 calculados mediante la estimación de máxima verosimilitud y el método de Newton-Raphson son los más adecuados, especialmente los obtenidos teniendo en cuenta los datos censurados. Tabla 4.12 Valores finales de τ0. m1 2345678910 τ0[h] 1 009.26 977.21 992.92 977.21 987.61 977.21 984.99 977.21 983.42 977.21 4.3 Caso de estudio 2: Fallos aleatorios Por definición, el caso de fallos aleatorios, al igual que el de fallos infantiles, no da lugar al cálculo, ni mucho menos al estudio de los intervalos entre mantenimientos preventivos. No obstante, se incluyen dos estudios con conjuntos de datos situados en cada una de esas regiones de la curva de bañera, con el fin de estudiar el desempeño de los métodos de estimación de los parámetros de Weibull. En particular, esta sección incluye el estudio del caso de fallos aleatorios. Como se puede observar en las tablas 4.13 y 4.14, los datos de fallo se encuentran repartidos a lo largo de la vida útil del equipo, registrándose fallos tanto al principio de su desempeño, como al final. Referente a los datos censurados, se observa que un equipo ha alcanzado τy se ha enviado al mantenimiento preventivo. Tabla 4.13 Datos de fallo para el caso de fallos aleatorios. Horas de funcionamiento antes del fallo 2 13 39 51 60 121 145 167 193 217 239 265 287 309 337 361 383 409 433 455 477 505 527 553 573 601 621 647 673 697 721 743 765 793 817 843 865 887 909 937 959 985 1 009 1 033 1 058 1 081 1 105 1 129 1 153 1 177 1 201 1 225 1 249 1 273 1 297 1 319 1 345 1 369 1 389 1 417 1 441 1 465 1 489 1 513 1 537 1 561 1 585 1 609 1 633 1 657 1 681 1 703 1 727 1 754 1 756 1 771 1 781 1 831 1 881 1 916 1 956 1 992 1999 Tabla 4.14 Datos censurados para el caso de fallos aleatorios. Horas de funcionamiento antes del fin de la observación 20 81 124 176 228 280 332 384 436 488 540 592 644 696 748 800 852 904 956 1 008 1 060 1 112 1 164 1 216 1 268 1 320 1 372 1 424 1 476 1 528 1 580 1 632 1 684 1 736 1 788 1 840 1 892 2 000 4.3 Caso de estudio 2: Fallos aleatorios 35 Al analizar este conjunto de datos con el método la regresión de rango de medianas, se obtienen los resultados plasmados en la figura 4.9. Figura 4.9 Resultados obtenidos mediante MRR para el caso de fallos aleatorios. Como se puede observar, se ha obtenido un valor de α de 1.01, muy cercano a la unidad. En este caso, la parábola posee un término cuadrático positivo, por lo que el valor de γ que lo anula debe de ser negativo, lo que va en contra de las restricciones impuestas en este proyecto, ya que γ>0 . Por ello, no tiene sentido calcular γpara el caso de fallos aleatorios. Respecto a la estimación de máxima verosimilitud, se ha añadido, en este caso, la restricción γ=0 en el complemento Solver, debido a la conclusión alcanzada sobre la vida garantizada al analizar los resultados de la regresión del rango de medianas. Figura 4.10 Iteraciones del método de Newton para el caso de fallos aleatorios. Por su parte, no puede llevarse a cabo el estudio de convergencia del algoritmo de Newton-Raphson al suponerse nulo el valor de la vida garantizada. De cualquier modo, se incluye la figura 4.10 para ilustrar el proceso iterativo que ha llevado a cabo el algoritmo. Como se puede observar en ella, el método converge en la iteración número 21. Por otro lado, de los estudios de la sección 4.2 se sacó como conclusión que existe cierta diferencia entre los resultados de cada método debido a que cada uno trata los datos de forma diferente y que algunos son más precisos que otros. En la tabla 4.15 se agrupan los diferentes resultados obtenidos por cada método en 36 Capítulo 4. Aplicación matemática este estudio. Tabla 4.15 Resumen de los resultados obtenidos para el caso de fallos aleatorios. Método α β γ MLE 1.51 1 376.15 0 Newton-Raphson 1.51 1 376.15 0 MRR 1.01 1 746.06 0 Al igual que en el análisis de fallo por desgaste, se registran valores más altos de α y más bajos de β en los métodos que emplean la función de verosimilitud. No obstante, en esta situación la diferencia es más acusada, de forma que el valor de α según la estimación de máxima verosimilitud y el método de Newton superan el valor unidad, entrando en la fase de desgaste. Esta mayor diferencia de cálculo puede deberse, a parte de a las diferencias relacionadas con el método de cálculo de cada método, al caso que se está estudiando. El caso de fallos aleatorios no suele darse en la vida real, donde predomina el desgaste. Por ello, al ser un caso especial, los métodos se muestran más sensibles, dando lugar a mayores diferencias. Por último, queda comentar que, al ser la regresión del rango de medianas un método más manejable, es más sencillo lograr ajustar una muestra de datos aleatoria con ella. La figura 4.11 incluye los valores de τ0 que se obtendrían con la expresión 3.17 al introducir los valores obtenidos con la regresión del rango de medianas, ya que es el método con el que se ha logrado alcanzar un α≈1. Figura 4.11 Intervalo óptimo entre mantenimientos preventivos según MRR para el caso de fallos aleatorios. Como se puede observar, la expresión 3.17 devuelve unos resultados muy grandes, haciendo imposible hablar de intervalos entre mantenimientos preventivos, lo cual es coherente con la situación en la que se encuentra la muestra analizada. 4.4 Caso de estudio 3: Fallos infantiles Finalmente, se presenta un conjunto de datos donde los fallos se dan a las pocas horas de comenzar a operar. Esto se puede ver reflejado en la tabla 4.16, donde más de la mitad de los fallos se producen antes de las 100 horas de funcionamiento. Del mismo modo, la tabla 4.17 agrupa al conjunto de datos censurados de este estudio. En este caso, un equipo logró alcanzar el mantenimiento preventivo. 4.4 Caso de estudio 3: Fallos infantiles 37 Tabla 4.16 Datos de fallo para el caso de fallos infantiles. Horas de funcionamiento antes del fallo 30 30 32 34 36 42 42 43 45 48 49 49 51 51 51 52 56 57 60 61 61 63 63 64 65 68 69 70 71 72 72 74 74 75 78 82 82 85 88 90 94 97 99 100 102 104 104 107 108 109 110 115 117 119 121 125 127 127 132 132 139 140 140 142 143 146 146 149 151 152 279 406 490 502 621 752 891 1 025 1 321 1 430 1 490 1 590 1 702 Tabla 4.17 Datos censurados para el caso de fallos infantiles. Horas de funcionamiento antes del fin de la observación 80 90 93 109 131 134 135 140 143 150 156 170 178 188 198 201 208 212 218 500 512 649 842 900 1 021 1 054 1 056 1 125 1 143 1 256 1 271 1 321 1 562 1 632 1 689 1 703 1 900 2 000 Según lo indicado en la sección 3.3, los conjuntos de datos con predominio de fallos prematuros deben de registrar un valor de α inferior a la unidad. Para comprobarlo, se hace uso de los mismos tres métodos que se han empleado a lo largo de este proyecto. Por su parte, los resultados obtenidos por medio de la regresión del rango de medianas se ven representados en la figura 4.12. Figura 4.12 Resultados obtenidos mediante MRR para el caso de fallos infantiles. Como se puede observar, se obtiene un valor de α de 0.92 y un valor de β de 378.87. En este caso, el valor de γ que anula el término cuadrático debe de ser muy pequeño para que no se produzcan errores en los cálculos. Como no se puede encontrar dicho valor, se asume nulo, siendo necesario volver a imponer, en este caso de estudio, la restricción γ=0en el complemento Solver de Excel. De nuevo, es imposible realizar el estudio de velocidad de convergencia de Newton-Raphson, convergiendo el algoritmo en la décimo tercera iteración. 44 Bibliografía [15] Didier Chauveau y Jean Diebolt Gilles Celeux, Stochastic versions of the EM algorithm: an experimental study in the mixture case, Journal of Statistical Computation and Simulation 55 (1996), no. 4, 287–314. [16] A. K. Gupta, Estimation of the mean and standard deviation of a normal population from a censored sample, Biometrika 39 (1952), no. 3/4, 260–273. [17] Hideo Hirose, Bias correction for the maximum likelihood estimates in the two-parameter Weibull distribution, IEEE Transactions on Dielectrics and Electrical Insulation 6(1999), 66–68. [18] J. R. M. Hosking, L-moments: Analysis and estimation of distributions using linear combinations of order statistics, Journal of the Royal Statistical Society. Series B (Methodological) 52 (1990), no. 1, 105–124. [19] James K. Starling, Chrsitina Mastrangelo y Youngjun Choe, Improving Weibull distribution estimation for generalized type I censored data using modified SMOTE, Reability Engineering and System Safety 221 (2021), 107505. [20] John P. Klein y Melvin L. Moeschberger, Survival analysis: Techniques for censored and truncated data, segunda ed., Springer-Verlag New York, Inc, 2003. [21] Kailash C. Kapur y Leonard R. Lamberson, Reliability in engineering design, primera ed., John Wiley & Sons, Inc, 1991. [22] M. Xie y L. C. Tang L. F. Zhang, Bias correction for the least squares estimator of Weibull shape parameter with complete and censored data, Reliability Engineering System Safety - RELIAB ENG SYST SAFETY 91 (2006), 930–939. [23] Kenneth L. Lange, A gradient algorithm locally equivalent to the EM algorithm, Journal of the Royal Statistical Society: Series B (Methodological) 57 (1995), 425–437. [24] Jerald F. Lawless, Statical models and methods for lifetime data, segunda ed., John Wiley & Sons, Inc, 2003. [25] Maha Alqam, Richard M. Bennett y Abdul-Hamid Zureick, Three-parameter vs. two-parameter Weibull distribution for pultruded composite material properties, Composite Structures 58 (2002), no. 4, 497– 503. [26] Mahdi Teimouri, Seyed M. Hoseini y Saralees Nadarajah, Comparison of estimation methods for the Weibull distribution, Statistics 47 (2013), no. 1, 93–109. [27] Barry Kurt Moser, Linear models: A mean model approach (a volume in probability and mathematical statistics), primera ed., Academic Press, 1996. [28] David M. Nanang, Suitability of the normal, log-normal and Weibull distributions for fitting diameter distributions of neem plantations in northern ghana, Forest Ecology and Management 103 (1998), 1–7. [29] Wayne Nelson, Applied life data analysis, primera ed., John Wiley & Sons, Inc, 1982. [30] Søren Feodor Nielsen, The stochastic EM algorithm: Estimation and asymptotic results, Bernoulli 6 (2000), no. 3, 457–489. [31] Milton Ohring, Engineering materials science, primera ed., Academic Press, noviembre 1995. [32] Philip Kosky, Robert Balmer, William Keat y George Wise, Exploring engineering. an introduction to engineering and design, quinta ed., Academic Press, abril 2020. [33] R. B. Abernethy, J. E. Breneman, C. H. Medlin y G. L. Reinman, Weybull analysis handbook, Tech. report, Defense Technical Information Center, noviembre 1933. [34] Rakesh Ranjan, Rijji Sen y Satyanshu K. Upadhyay, Bayes analysis of some important lifetime models using MCMC based approaches when the observations are left truncated and right censored, Reliability Engineering and System Safety 214 (2021), 107747. Bibliografía 45 [35] Ricardo Medel Esquivel, Isidro Gomez-Vargas, J. Alberto Vazquez y Ricardo García-Salcedo, An introduction to Markov chain Monte Carlo, Boletín de Estadistica e Investigacion Operativa 37 (2021), 47–84. [36] Riffat Jabeen, Azaz Ahmad, Navid Feroze y Ghausia M. Gilani, Estimation of location and scale parameters of Weibull distribution using generalized order statistics under type II singly and doubly censored data, International Journal of Advanced Science and Technology 55 (2013), 67–80. [37] Rola M. Musleh y Amal Helu, Estimation of the inverse Weibull distribution based on progressively censored data: Comparative study, Reliability Engineering & System Safety 131 (2014), 216–227. [38] Robert Ross, Formulas to describe the bias and standard deviation of the ML-estimated Weibull shape parameter, IEEE Transactions on Dielectrics and Electrical Insulation 1(1994), 247–253. [39] ,Bias and standard deviation due to Weibull parameter estimation for small data sets, IEEE Transactions on Dielectrics and Electrical Insulation 3(1996), 28–42. [40] Russell B. Millar, Maximum likelihood estimation and inference: with examples in R, SAS and ADMB, primera ed., John Wiley & Sons, Inc, 2011. [41] S. B. Kang, Y. S. Cho y S. H. Choi, Approximate maximum likelihood estimation for the three-parameter Weibull Distribution, The Korean Communications in Statistics 8(2001), no. 1, 209–217. [42] Stephen Shifley and Ellen Lentz, Quick estimation of the three-parameter Weibull to describe tree size distributions, Forest Ecology and Management 13 (1985), 195–203. [43] David Smith, Reliability, maintainability and risk. practical methods for engineers including reliability centred maintenance and safety-related systems, séptima ed., Butterworth-Heinemann, abril 2005. [44] Steven Brown, Beata Walczak y Romà Tauler, Comprehensive chemometrics: Chemical and biochemical data analysis, primera ed., Elsevier Science, Ltd, 2009. [45] Jeannet Gunneklev Støtvig, Censored Weibull distributed data in experimental design, Ph.D. thesis, Norwegian University of Science and Technology, 2014. [46] Switamy Angnitha Purba, Sutarman y Open D. Sembiring , Maximum likelihood based on Newton- Raphson, Fisher scoring and expectation maximization algorithm application on accident data., International Journal of Advanced Research 6(2018), 965–969. [47] William Q. Meeker y Luis A. Escobar, Statical methods for reliability data, primera ed., John Wiley & Sons, Inc, 1998. [48] N. Balakrishnan y A. Clifford Cohen, Order statistics and inference: Estimation methods, primera ed., Academic Press, INC., 1991. [49] Patrick D. T. O’Connor y André Kleyner, Practical reliability engineering, quinta ed., John Wiley & Sons, Inc, 2012. [50] Kandethody M. Ramachandran y Chris P. Tsokos, Mathematical statistics with applications in R, tercera ed., Elsevier Inc, 2020. [51] N. Balakrishnan y Debanjan Mitra, Left truncated and right censored Weibull data and likelihood inference with an illustration, Computational Statistics and Data Analysis 56 (2012), no. 12, 4011–4025. [52] Mostafa Bacha y Gilles Celeux, Bayesian estimation of a Weibull distribution in a highly censored and small sample setting, Research Report RR-2993, INRIA, 1996, Projet ISIS. [53] Jan Merganič y Hubert Sterba, Characterisation of diameter distribution using the Weibull function: Method of moments, European Journal of Forest Research 125 (2006), 427–439. [54] G. Celeux y J. Diebolt, The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem, Computational Statistics Quarterly 2(1986), 73–82. 46 Bibliografía [55] N. Balakrishnan y J. Varadan, Approximate MLEs for the location and scale parameters of the extreme value distribution with censoring, IEEE Transactions on Reliability 40 (1991), no. 2, 146–151. [56] Luis Ferreira y José Silva, Parameter estimation for Weibull distribution with right censored data using EM algorithm, Eksploatacja i Niezawodnosc - Maintenance and Reliability 19 (2017), 310–315. [57] Florence Ducros y Patrick Pamphile, Bayesian estimation of Weibull mixture in heavily censored data setting, Reliability Engineering and System Safety 180 (2018), no. C, 453–462. [58] Takeshi Emura y Shau-Kai Shiu, Estimation and model selection for left-truncated and right-censored lifetime data with application to electric power transformers analysis, Communications in Statistics - Simulation and Computation 45 (2016), no. 9, 3171–3189. [59] Stanley J. Zarnoch y Tommy R. Dell, An evaluation of percentile and maximum likelihood estimators of weibull paremeters, Forest Science 31 (1985), 260–268. [60] Ulrike Genschel y William Meeker, A comparison of maximum likelihood and median-rank regression for Weibull estimation, Quality Engineering 22 (2010), 236–255. [61] Anwar Hossain y William Zimmer, Comparison of estimation methods for Weibull parameters: Complete and censored samples, Journal of Statistical Computation and Simulation 73 (2003), no. 2, 145–153. [62] Stelios H. Zanakis, A simulation study of some simple estimators for the three-parameter Weibull distribution, Journal of Statistical Computation and Simulation 9(1979), 101–116.