1.4 Modelos Probabilísticos em Confiabilidade
No tópico anterior apresentamos um estudo envolvendo um teste de vida com as válvulas do Exemplo 2.1. O fabricante estava interessado em estimar as seguintes características do seu produto:
O tempo médio de falha (MTTF);
O tempo no qual 10% dos mecanismos estarão fora de operação;
O percentual esperado de falhas nos dois primeiros anos de uso.
As duas últimas características foram estimadas usando os estimadores tabela de vida e Kaplan-Meier. Essas técnicas são chamadas de não paramétricas, pois não necessitam da especificação de nenhuma distribuição de probabilidade para o tempo de vida do componente. Entretanto, com essas técnicas muitas vezes não é possível obter uma boa estimativa de um parâmetro importante: o tempo médio de falha, MTTF. No entanto, existem técnicas estatísticas disponíveis que requerem a especificação de uma distribuição de probabilidade e que torna a estimação desse parâmetro (MTTF) mais confiável. Estas técnicas são chamadas de paramétricas. Quando a distribuição de probabilidade escolhida está adequada aos dados, as estimativas paramétricas são mais confiáveis que as estimativas não paramétricas.
O objetivo deste tópico é apresentar métodos estatísticos paramétricos para a análise do tempo de vida, baseados em modelos probabilísticos importantes, tais como o Exponencial, Weibull, Gumbel ou Valor Extremo e Log-normal, que são discutidos a seguir.
As duas últimas características foram estimadas usando os estimadores tabela de vida e Kaplan-Meier. Essas técnicas são chamadas de não paramétricas, pois não necessitam da especificação de nenhuma distribuição de probabilidade para o tempo de vida do componente. Entretanto, com essas técnicas muitas vezes não é possível obter uma boa estimativa de um parâmetro importante: o tempo médio de falha, MTTF. No entanto, existem técnicas estatísticas disponíveis que requerem a especificação de uma distribuição de probabilidade e que torna a estimação desse parâmetro (MTTF) mais confiável. Estas técnicas são chamadas de paramétricas. Quando a distribuição de probabilidade escolhida está adequada aos dados, as estimativas paramétricas são mais confiáveis que as estimativas não paramétricas.
O objetivo deste tópico é apresentar métodos estatísticos paramétricos para a análise do tempo de vida, baseados em modelos probabilísticos importantes, tais como o Exponencial, Weibull, Gumbel ou Valor Extremo e Log-normal, que são discutidos a seguir.
4.1 - Modelos Probabilísticos para o Tempo de Falha
Existe uma série de modelos probabilísticos utilizados na análise de dados de confiabilidade, alguns destes modelos ocupam uma posição de destaque por sua comprovada adequação a várias situações práticas. Entre esses modelos podemos citar o Exponencial, Weibull, Valor Extremo ou Gumbel e o Log-normal.
Figura 1.4.1: Principais modelos utilizados em análise de confiabilidade
É importante entender que cada distribuição de probabilidade pode gerar diferentes estimativas para as características de durabilidade do produto. Dessa forma, a utilização de um modelo inadequado levará a erros grosseiros de estimativas. Portanto, a escolha de um modelo adequado para descrever o tempo de falha de um determinado produto deve ser feita com bastante cuidado.
4.1.1 - Distribuição Exponencial
A distribuição exponencial se caracteriza por uma função taxa de falha constante, sendo a única distribuição absolutamente contínua com essa propriedade. Ela é considerada uma das mais simples em termos matemáticos. Esta distribuição tem sido extensivamente utilizada para modelar o tempo de vida de certos produtos e materiais, tais como óleos isolantes, dielétricos, entre outros. A função densidade de probabilidade para um tempo de falha T com distribuição exponencial é dada por
$$f(t)= \alpha e^{-\alpha t}, ~~~~~t~>~0,$$
sendo $1/\alpha~>~0$ o tempo médio de vida (MTTF), por isso o parâmetro $\alpha$ tem a mesma unidade do tempo de vida. Isto é, se o tempo é medido em horas, o valor de $1/\alpha$ representa o tempo médio em horas. Algumas formas típicas da função densidade de probabilidade da distribuição exponencial são mostradas na Figura 4.1.1.1.
Figura 1.4.2: Funções densidade de probabilidade da distribuição exponencial.
A função de confiabilidade R(t), que é a probabilidade do produto continuar funcionando além do tempo $t$ é dada por
$$R(t)= 1- F(t) = e^{-\alpha t},~~~~~t~>~0.$$
A Figura 4.1.1.2 mostra formas típicas da função confiabilidade da distribuição exponencial.
Figura 1.4.3: Funções de confiabilidade para a distribuição exponencial.
A função taxa de falha associada a distribuição Exponencial é constante e é dada por
$$h(t) = \alpha,~~~~~t~>~0,$$
isto é, tanto uma unidade que está em operação a 20 horas quanto uma unidade que está em operação a 40 horas tem a mesma chance de falharem em um intervalo futuro de mesmo comprimento. Esta propriedade é chamada de falta de memória da distribuição.
Outras características importantes que dizem respeito a durabilidade de um componente são: a média (MTTF), a variância e os percentis da distribuição do tempo de falha. A média da distribuição exponencial (MTTF) e a variância são dadas por
$$E(T)= \hbox{MTTF}=\dfrac{1}{\alpha} \qquad \hbox{e} \qquad Var(T)=\dfrac{1}{\alpha^2}.$$
O quantil 100×p% corresponde ao tempo esperado em que 100×p% dos produtos falham, sendo obtido a partir da função de confiabilidade como
$$t_{p} = - \dfrac{1}{\alpha} \log(1-p) \tag{4.1.1.1}$$
Muitas vezes, em estudos de durabilidade queremos conhecer percentis pequenos como 1%, que informam sobre tempos de falhas prematuras. Outro quantil importante é a mediana ou quantil 50%, que informa sobre o tempo em que metade das unidades falham. A média da distribuição exponencial corresponde ao quantil 63%.
Exemplo 4.1.1.1
Considere que o tempo até a falha do ventilador de motores a diesel segue uma distribuição exponencial com $\alpha=1/28.700$ horas. Vamos calcular a probabilidade de um destes ventiladores não falhar nas primeiras 8.000 horas de funcionamento. Dessa forma, usando a função de confiabilidade R(t) para t = 8.000 temos
$$R(8.000) = \exp\left(-\dfrac{8.000}{28.700}\right) = 0,76.$$
Se 8.000 horas for o tempo de garantia dado pelo fabricante, então espera-se que 24% dos ventiladores falhem antes do término da garantia.
O quantil 1% é obtido a partir da equação (4.1.1.1) como
$$t_{0,01} = -28.700 \times \log\left(1- 0,01\right) = 288,44.$$
Logo, espera-se que 1% dos ventiladores falhem nas primeiras 288 horas de uso. De forma similar, obtemos a mediana como sendo 19.893 horas.
Para uma visão geral das características da distribuição exponencial utilizadas em confiabilidade você pode utilizar também o Software Action, mais especificamente a ferramenta Overview para confiabilidade.
4.1.2 - Distribuição de Weibull
A distribuição de Weibull foi proposta originalmente por W. Weibull em 1954 em estudos relacionados ao tempo de falha devido a fadiga de metais. Ela é frequentemente usada para descrever o tempo de vida de produtos industriais. A sua popularidade em aplicações práticas deve-se ao fato dela apresentar uma grande variedade de formas, todas com uma propriedade básica: a sua função de taxa de falha é monótona, isto é, ela é estritamente crescente, estritamente decrescente ou constante. Ela descreve adequadamente a vida de mancais, componentes eletrônicos, cerâmicas, capacitores e dielétricos. Sua função densidade de probabilidade é dada por
$$f(t) = \dfrac{\delta}{\alpha^\delta}~t^{\delta-1}\exp \left[-\left(\dfrac{t}{\alpha}\right)^{\delta} \right]~~~~t~>~0~, \tag{4.1.2.1}$$
sendo $\delta~>~0$ e $\alpha~>~0$ os parâmetros de forma e de escala, respectivamente.
Observe que se tomarmos o parâmetro de forma $\delta = 1$ na equação (4.1.2.1), obtemos a densidade de probabilidade da Exponencial com parâmetro $\alpha$, isto é, a distribuição de Weibull é uma generalização da distribuição Exponencial.
Figura 1.4.4: Funções densidade de probabilidade da distribuição de Weibull com $\alpha = 1.$
A função de confiabilidade da distribuição de Weibull é dada por
$$R(t) = \exp\left(-\left(\dfrac{t}{\alpha}\right)^{\delta}\right)~~~~t~>~0~ \tag{4.1.2.2}$$
Figura 1.4.5: Funções de confiabilidade para a distribuição de Weibull com $\alpha=1.$
A função de risco é dada por
$$h(t) = \dfrac{\delta}{\alpha}\left(\dfrac{t}{\alpha}\right)^{\delta-1},~~~~~t~>~0.$$
A Figura 4.1.2.3 mostra algumas formas da função de risco para a distribuição de Weibull. Observe que h(t) é estritamente crescente para $\delta~>~1$ e estritamente decrescente para $\delta~<~1$. Como a distribuição exponencial é um caso particular da distribuição de Weibull quando $\delta = 1$, a taxa de falha fica constante neste caso.
Figura 1.4.6: Taxa de falha da distribuição de Weibull com $\alpha=1.$
O tempo médio de vida (MTTF) e a variância da distribuição de Weibull são obtidos pelas equações
$$\hbox{MTTF}=\alpha~\Gamma\left[1 + \dfrac{1}{\delta}\right] \tag{4.1.2.3}$$
$$\hbox{Var}(T)=\alpha^2~\left(\Gamma\left[1+\dfrac{2}{\delta}\right]-\Gamma^2\left[1+\dfrac{1}{\delta}\right]\right),$$
sendo $\Gamma(\cdot)$ denominada função Gamma e definida como
$$\Gamma(x) = \int_0^{\infty} u^{x-1}~e^{-u}du,~~~~~~~x \geq 0$$
Quando x é um número natural positivo, temos que
$$\Gamma(x) = (x - 1)!$$
O quantil 100×p% é obtido a partir da equação (4.1.2.2) como
$$t_p = \alpha\left[-\log(1-p)\right]^{1/\delta} \tag{4.1.2.4}.$$
Exemplo 4.1.2.1
Suponha que o tempo de vida de um capacitor obedeça a uma distribuição de Weibull com parâmetros $\alpha = 100.000$ e $\delta = 0,5$.
A confiabilidade para um ano, t = 8.760, é obtida a partir da equação (4.1.2.2) como
$$R(t)=\exp\left[-\left(\dfrac{8.760}{100.000}\right)^{0,5}\right] = \exp\left[-0,2960\right] = 0,7438,$$
isto é, a probabilidade do capacitor operar por mais de um ano é de 74%.
O tempo médio de vida (MTTF) deste capacitor é obtido a partir da equação (4.1.2.3) como
$$\hbox{MTTF}=100.000 \ast \Gamma\left(1+\dfrac{1}{0,5}\right)=100.000 \ast \Gamma(3) = 200.000,$$
isto é, em média esse capacitor dura 200.000 horas.
Da equação (4.1.2.4) obtemos o quantil 10% como
$$t_{0,10}=100.000 \ast \left[-\log(0,9)\right]^{1/0,5}=100.000 \ast \left[0,1054\right]^{2} = 1.110,92,$$
isto é, espera-se que 10% das unidades falhem antes de atingir 1.111 horas de uso.
Para uma visão geral das características da distribuição Weibull em confiabilidade você pode utilizar também o Software Action, mais especificamente a ferramenta Overview para confiabilidade.
4.1.3 - Distribuição de Valor Extremo (ou de Gumbel)
Neste tópico, introduzimos uma distribuição que é bastante relacionada à distribuição de Weibull. Ela é chamada de distribuição de Valor Extremo (ou de Gumbel) e surge quando se considera o logaritmo natural de uma variável aleatória com distribuição de Weibull. Isto é, se uma variável T tem distribuição de Weibull com parâmetros $\alpha$ (escala) e $\delta$ (forma), então a variável Y = log(T) tem uma distribuição de Valor Extremo com função densidade de probabilidade dada por
$$f(y)= \dfrac{1}{\sigma}\exp\left[\dfrac{y-\mu}{\sigma}-\exp\left(\dfrac{y-\mu}{\sigma}\right)\right],~~~~~~\infty~<~y~<~\infty \tag{4.1.3.1}$$
sendo $\mu = \log(\alpha )$ e $ \sigma = \dfrac{1}{\delta}$ os parâmetros de locação e escala, respectivamente.
Se $\mu = 0$ e $\sigma = 1$ temos a distribuição do Valor Extremo padrão.
Figura 1.4.7: Funções densidade de probabilidade da distribuição de Valor Extremo com $\mu=0.$
A função de confiabilidade da variável Y é dada por
$$R(y) = \exp \left[-\exp \left(\dfrac{y-\mu}{\sigma}\right)\right], \qquad -\infty~<~y~<~\infty \tag{4.1.3.2}$$
Figura 1.4.8: Funções de confiabilidade para a distribuição de Valor Extremo para $\mu=0.$
Considerando as funções dadas pelas equações (4.1.3.1) e (4.1.3.2), segue que a função de risco da distribuição de Valor Extremo é dada por
$$h(y) = \dfrac{1}{\sigma} \exp \left[ \dfrac{y - \mu}{\sigma}\right], \qquad -\infty~<~y~<~\infty \tag{4.1.3.3}.$$
Figura 1.4.9: Taxas de falha da distribuição de Valor Extremo com $\mu=0.$
O tempo médio de vida (MTTF) e a variância de Y são dados por
$$\hbox{MTTF} = \mu - \gamma , \sigma$$
$$\hbox{Var}(Y) = \dfrac{\pi^2 , \sigma^2}{6},$$
sendo $\gamma \simeq 0,5772$ denominada constante de Euler.
O quantil 100×p% é dado por
$$t_p = \mu + \sigma \log \left[-\log(1-p) \right].$$
Na análise de dados de durabilidade muitas vezes é conveniente trabalhar com o logaritmo dos tempos de vida observados. Esse fato será bastante explorado na Seção 6, em que os modelos de regressão são usados para analisar os dados provenientes de testes de vida acelerados. Dessa forma, se os dados tiverem uma distribuição de Weibull, a distribuição de Valor Extremo aparecerá naturalmente na modelagem.
Para uma visão geral das características da distribuição de valor extremo utilizadas em confiabilidade você pode utilizar também o Software Action, mais especificamente a ferramenta Overview para confiabilidade.
4.1.4 - Distribuição Log-Normal
Assim como a distribuição de Weibull, a distribuição log-normal é muito usada para caracterizar tempo de vida de produtos e materiais. Isso inclui fadiga de metal, semicondutores, diodos e isolação elétrica.
A função densidade de probabilidade da distribuição log-normal é dada por
$$f(t) = \dfrac{1}{\sqrt{2 \pi}, , t \sigma} , \exp \left[-\dfrac{1}{2} \left(\dfrac{\log (t) - \mu}{\sigma} \right)^2\right] ,\qquad t~>~0, \tag{4.1.4.1}$$
sendo $-\infty~<~\mu~<~\infty$ e $\sigma~>~0$.
Figura 1.4.10: Funções densidade de probabilidade da distribuição log-Normal com $\mu=0.$
Existe uma relação entre as distribuições log-normal e normal, similar à existente entre as distribuições de Weibull e de valor extremo. O logaritmo de uma variável que segue distribuição log-normal com parâmetros $\mu$ e $\sigma$ tem distribuição normal com média $\mu$ e desvio-padrão $\sigma$. Essa relação significa que dados provenientes de uma distribuição log-normal podem ser analisados segundo uma distribuição normal, se considerarmos o logaritmo dos dados ao invés dos valores originais.
A função de confiabilidade de uma distribuição log-normal é dada por
$$R(t) = \Phi \left(-\dfrac{\log(t)-\mu}{\sigma} \right), \tag{4.1.4.2}$$
em que $\Phi(.)$ é a função de distribuição acumulada de uma distribuição normal padrão.
Figura 1.4.11: Funções de confiabilidade para a distribuição log-Normal com $\mu=1.$
A taxa de falha é dada por
$$h(t) = \dfrac{f(t)}{R(t)}, \qquad t~>~0 \tag{4.1.4.3}$$
Figura 1.4.12: Taxas de falha da distribuição log-Normal com $\mu=1.$
O tempo médio de vida e a variância da distribuição log-normal são dados, respectivamente, por
$$\hbox{MTTF} = \exp \left(\mu + \dfrac{\sigma^2}{2} \right), \tag{4.1.4.4}$$
$$Var(T) = \exp (2 \mu + \sigma^2) [\exp(\sigma^2) - 1]. \tag{4.1.4.5}$$
O quantil 100×p% da distribuição Log-normal é dado pela expressão
$$t_p = \exp(z_p \sigma + \mu), \tag{4.1.4.6}$$
sendo $z_p$ o quantil 100×p% da distribuição normal padrão.
Exemplo 4.1.4.1
Para ilustrar o uso da distribuição log-normal, considere o tempo de vida de isolações da classe H. Supomos que o tempo de vida de uma isolação, na temperatura de uso, tem uma distribuição log-normal com parâmetros $\mu = 9,65$ e $\sigma = 0,1053$.
A confiabilidade de uma isolação nas 20.000 primeiras horas de uso é obtida a partir da equação (4.1.4.2) como
$$R(20.000) = \Phi \left( -\dfrac{\log(20.000) - 9,65}{0,1053}\right) = \Phi (-2,4073) = 0,008.$$
Isso significa que a grande maioria (99,2%) das isolações falham nas 20.000 primeiras horas de uso.
O tempo médio de vida (MTTF) de uma isolação é obtido a partir da equação (4.1.4.4) como
$$\hbox{MTTF} = \exp \left(9,65 + \dfrac{0,1053^2}{2}\right) = \exp \left(9,6555\right) = 15.607,39.$$
A partir da equação (4.1.4.6) podemos obter a mediana, ou seja, o quantil 50% da distribuição como sendo
$$t_{0,5} = \exp (0 \times 0,1053 + 9,65) = \exp (9,65) = 15.521,79.$$
Para uma visão geral das características da distribuição log-normal utilizadas em confiabilidade você pode utilizar também o Software Action, mais especificamente a ferramenta Overview para confiabilidade.
4.2 - Estimando os Parâmetros dos Modelos
Os modelos probabilísticos apresentados na seção anterior são caracterizados por quantidades desconhecidas, denominadas parâmetros. Por exemplo, os modelos de Weibull e Log-normal são caracterizados por dois parâmetros enquanto que o modelo Exponencial é caracterizado por apenas um. Essas quantidades conferem uma forma geral aos modelos probabilísticos. Entretanto, em cada estudo de confiabilidade, os parâmetros devem ser estimados a partir das observações amostrais para que o modelo fique determinado e assim possa ser possível responder às perguntas de interesse. Existem alguns métodos de estimação conhecidos na literatura estatística. O mais conhecido é o método de mínimos quadrados, geralmente apresentado em cursos básicos de estatística dentro do contexto de regressão linear. No entanto, esse método é inapropriado para estudos de durabilidade. A principal razão é a sua incapacidade de incorporar censuras no seu processo de estimação. Sendo assim, o método de máxima verossimilhança surge como uma opção apropriada para esse tipo de dados. Esse método incorpora censuras, sendo relativamente simples de ser entendido e com propriedades ótimas para grandes amostras.
4.2.1 - Método de Máxima Verossimilhança
O método de máxima verossimilhança trata o problema de estimação da seguinte forma: baseado nos resultados obtidos pela amostra, devemos determinar qual a distribuição, dentre todas aquelas definidas pelos possíveis valores de seus parâmetros, com maior possibilidade de ter gerado tal amostra. Em outras palavras, se por exemplo a distribuição do tempo de falha é a de Weibull, para cada combinação diferente de $\alpha$ e $\delta$ tem-se diferentes distribuições de Weibull. O estimador de máxima verossimilhança escolhe aquele par de $\alpha$ e $\delta$ que melhor explica a amostra observada. A seguir discutiremos as idéias do método de máxima verossimilhança para conceitos matemáticos a partir dos quais será possível obter estimadores para os parâmetros. Suponha uma amostra de observações $t_{1}$, …, $t_{n}$ de uma certa população de interesse. Considere inicialmente que todas as observações são não-censuradas. A população é caracterizada pela sua função de densidade f(t). Por exemplo, se $f(t) = \alpha \exp(-\alpha t)$, significa que as observações vem de uma distribuição Exponencial com um parâmetro a ser estimado. A função de verossimilhança para um parâmetro genérico é dada por
$$L(\theta)= \displaystyle\prod_{i=1}^{n}f(t_i;\theta) \tag{4.2.1.1}$$
Note na expressão (4.2.1.1) que θ pode representar um único parâmetro ou um conjunto de parâmetros. Por exemplo, no modelo Log-normal temos θ = (μ, σ). A tradução em termos matemáticos para a frase “a distribuição que melhor explica a amostra observada” é achar o valor θ que maximiza a função L(θ). Isto é, achar o valor de θ que maximiza a probabilidade da amostra observada ocorrer. A função de verossimilhança L(θ) mostra que a contribuição de cada observação não-censurada é a sua função de densidade. Por outro lado, a contribuição de cada observação censurada não é a sua função de densidade. Essas observações somente nos informam que o tempo de falha é maior do que o tempo de censura observado e portanto, sua contribuição para L(θ) é a sua função de confiabilidade R(t). As observações podem então ser divididas em dois conjuntos, as r primeiras são as não-censuradas (1, 2, …, r), e as n - r seguintes são as censuradas (r + 1, r + 2, …, n). Com isso, a função de verossimilhança assume a seguinte forma
$$L(\theta)=\prod_{i=1}^{r}f(t_i;\theta)\prod_{i=r+1}^{n}R(t_i;\theta) \tag{4.2.1.2}$$
A expressão (4.2.1.2) para a verossimilhança é válida para os mecanismos de censura do tipo I e II sob a suposição de que o mecanismo de censura é não-informativo, ou seja, não carrega informações sobre os parâmetros, a mesma vale também para o mecanismo do tipo aleatório. Essa suposição é razoável em estudos de durabilidade e é sempre conveniente trabalhar com o logaritmo da função de verossimilhança. Os estimadores de máxima verossimilhança são os valores θ que maximizam L(θ) ou equivalente log(L(θ)). Eles são encontrados resolvendo o sistema de equações
$$U(\theta)= \dfrac{\partial \log(L(\theta))}{\partial\theta}=0$$
Aplicações
Os cálculos realizados para obter os estimadores de máxima verossimilhança são ilustrados abaixo para as distribuições Exponencial e de Weibull. No caso da distribuição de Weibull, como não existem expressões fechadas para os estimadores de $\alpha$ e $\delta$ optamos por apresentar os passos seguidos pelo método numérico. Nesse caso, as estimativas para um conjunto de dados de durabilidade devem ser obtidas por meio de um pacote estatístico.
Daqui em diante, vamos supor que $t_{1}$, …, $t_{n}$ é uma amostra de observações independentes, em que $t_{1}$, …, $t_{r}$ são os tempos observados de falha e $t_{r+1}$, …, $t_{n}$ são observações censuradas.
Distribuição Exponencial
A função de verossimilhança para a distribuição Exponencial é obtida a partir da expressão (4.2.1.2) como
$$L(\alpha)=\left(\prod_{i=1}^{r}f(t_i;\theta)\right)\times\left(\prod_{i=r+1}^{n}R(t_i;\theta)\right)$$
$$=\left(\prod_{i=1}^{r}\alpha \exp(-\alpha t_{i})\right) \times\left(\prod_{i=r+1}^{n} \exp(-\alpha t_{i})\right)$$
$$=\alpha^{r} \exp\left(-\alpha \sum_{i=1}^{n} t_{i}\right).$$
Com isso, o logaritmo da função de verossimilhança é dado por
$$\log L(\alpha)=r \log (\alpha) - \alpha \displaystyle\sum_{i=1}^{n}t_{i}$$
Derivando essa expressão em relação a $\alpha$, obtemos
$$\dfrac{\partial \log L(\alpha)}{\partial \alpha}=\dfrac{r}{\alpha}-\sum_{i=1}^{n}t_{i},$$
e igualando a zero, temos que a expressão do estimador de máxima verossimilhança $\widehat{\alpha}$ é dada por
$$\widehat{\alpha}=\dfrac{r}{\displaystyle\sum_{i=1}^{n}t_i}.$$
O termo $\displaystyle \sum_{i=1}^{n} t_i$ é denominado “tempo total sob teste”. Observe que se todas as observações são não-censuradas, temos que a média amostral é dada por $\widehat{\alpha}=1/\overline{t}.$
Distribuição de Weibull
A função de verossimilhança de $\alpha$ e $\delta$ para a distribuição de Weibull é dada por
$$L(\alpha,\delta)=\left(\prod_{i=1}^{r}f(t_i;\theta)\right)\times\left(\prod_{i=r+1}^{n}R(t_i;\theta)\right)$$
$$=\left(\prod_{i=1}^{r} \dfrac{\delta}{\alpha^{\delta}} t_i^{\delta - 1} \exp \left[-\left(\dfrac{t_i}{\alpha}\right)^{\delta}\right]\right)\times\left(\prod_{i=r+1}^{n}\exp\left[-\left(\dfrac{t_i}{\alpha}\right)^{\delta}\right]\right)$$
$$=\dfrac{\delta^r}{\alpha^{r\delta}}\left(\prod_{i=1}^{r}t_i\right)^{\delta-1}\exp\left(-\dfrac{1}{\alpha^{\delta}} \sum_{i=1}^{n} t_i^{\delta}\right).$$
Com isso, a função de log-verossimilhança é dada por
$$\log L(\alpha,\delta)= r \log(\delta) - r \delta \log (\alpha) + (\delta-1)\sum_{i=1}^{r} \log(t_i) - \dfrac{1}{\alpha^{\delta}}\sum_{i=1}^{n} t_i^{\delta}$$
Derivando essa expressão em relação a $\alpha$ e $\delta$ e igualando a zero, obtemos as seguintes expressões para os estimadores de máxima verossimilhança $\widehat{\alpha}$ e $\widehat{\delta}$
$$\dfrac{\sum_{i=1}^{n}t_i^{\widehat{\delta}}\log(t_i)}{\sum_{i=1}^{n}t_i^{\widehat{\delta}}}- \dfrac{1}{\widehat{\delta}}-\dfrac{1}{r}\sum_{i=1}^{r}\log(t_i)=0, \tag{4.2.1.3}$$
$$\widehat{\alpha}=\left(\dfrac{1}{r}\sum_{i=1}^{n}t_i^{\widehat{\delta}}\right)^{1/\widehat{\delta}}. \tag{4.2.1.4}$$
Os estimadores de máxima verossimilhança são os valores $\widehat{\alpha}$ e $\widehat{\delta}$ que satisfazem as equações (4.2.1.3) e (4.2.1.4). A solução desse sistema de equações para um conjunto de dados particular deve ser obtida por meio de um método numérico. Aqui utilizaremos o método de Newton-Raphson que usa a matriz de derivadas segundas (F) da função de log-verossimilhança, sua expressão é dada por
$$\widehat{\theta}^{(k+1)}=\widehat{\theta}^{(k)}-F^{-1} (\widehat{\theta}^{(k)})U(\widehat{\theta}^{(k)}), \tag{4.2.1.5}$$
em que
$$U(\theta)=\dfrac{\partial \log L(\theta)}{\partial \theta}.$$
A expressão (4.2.1.5) é baseada na expansão de $U(\widehat{\theta}^{(k)})$ em série de Taylor em torno de $\widehat{\theta}^{(k)}$. Partindo de um valor inicial $\widehat{\theta}^{(0)}=0$ o método atualiza esse valor a cada passo, convergindo para a solução desejada. Em geral, obtemos convergência em poucos passos com um erro relativo, em média, menor que 0,001 entre dois passos consecutivos. Observe que a matriz de derivadas F para o modelo Exponencial se reduz a um único número, dado por
$$F(\alpha)=\dfrac{\partial^2 \log L(\alpha)}{\partial \alpha^2}=\dfrac{r}{\alpha^2}-2\dfrac{\sum_{i=1}^{n}t_i}{\alpha^3}.$$
Já para o modelo de Weibull $F(\alpha, \delta)$ é uma matriz simétrica 2x2 com os seguintes elementos
$$F_{11}(\alpha,\delta)=\dfrac{\partial^2 \log L(\alpha,\delta)}{\partial \alpha^2}$$
$$F_{22}(\alpha,\delta)=\dfrac{\partial^2 \log L(\alpha,\delta)}{\partial \delta^2}$$
$$F_{11}(\alpha,\delta)=F_{21}(\alpha,\delta)=\dfrac{\partial^2 \log L(\alpha,\delta)}{\partial \alpha \partial \delta}$$
Intervalos de Confiança para os parâmetros
O método de máxima verossimilhança foi usado para obter estimadores para os parâmetros do modelo. Esses valores são chamados de estimadores pontuais. Esse método também permite a construção de intervalos de confiança para os parâmetros. Isso é feito a partir das propriedades para grandes amostras desses estimadores. As justificativas matemáticas dessas propriedades são bastante complexas e nesse texto apresentaremos apenas as propriedades que são suficientes para os objetivos propostos. A propriedade ou resultado mais importante diz respeito à precisão do estimador de máxima verossimilhança e estabelece que
$$Var(\widehat{\theta}) \approx -[E(F(\theta))]^{-1}$$
Ou seja, isso significa que a matriz de variâncias-covariâncias do estimador de máxima verossimilhança é aproximadamente o negativo da inversa da matriz de segundas derivadas de $\log(L(θ))$ esperada. Nas situações em que a esperança é impossível ou difícil de ser calculada, usa-se simplesmente $-(F^{-1}(θ))$, em que $F(θ)$ é a matriz de segundas derivadas de $\log(L(θ))$. Os elementos da diagonal principal são as variâncias dos estimadores e os outros suas respectivas covariâncias. Geralmente $Var(θ)$ depende de $θ$. Uma estimativa para $Var(θ)$ é então obtida substituindo $θ$ por $\widehat{\theta}$. Na construção de intervalos de confiança é necessário ter uma estimativa para o erro-padrão de $\widehat{\theta}=\sqrt{Var(\widehat{\theta})}.$ No caso particular em que θ é um escalar, um intervalo de $(1-\alpha)100\char37$ de confiança para $θ$ é
$$\widehat{\theta} \pm Z_{\alpha/2}\sqrt{Var(\widehat{\theta})}.$$
Por exemplo, um intervalo de 95% de confiança para o parâmetro do modelo exponencial é dado por
$$\alpha \pm 1.96 \times \sqrt{\dfrac{\widehat{\alpha}^2}{r}}$$
pois,
$$E \left[\dfrac{r}{\alpha^2}-\dfrac{-2 \sum_{i=1}^{n}t_i}{\alpha^3} \right]=\dfrac{-r}{\alpha^2}$$
No caso em que $θ$ é um vetor de parâmetros, um intervalo de confiança pode ser construído para cada parâmetro separadamente. Basta obter uma estimativa para o erro-padrão a partir da matriz de variância-covariância $Var(\widehat{\theta}).$ Consideremos agora $\theta=(\delta, \alpha)$ como no modelo de Weibull. Algumas vezes o interesse pode ser estimar uma função dos parâmetros $\phi = g(\alpha ,\delta),$ por exemplo a mediana da distribuição de Weibull, $t_{0,5} = \alpha[-\log(1 - 0,5 )]^{1/\delta}.$ O estimador de máxima verossimilhança $\phi$ é $\widehat{\phi}=g(\widehat{\alpha},\widehat{\delta}),$ ou seja, para estimar $\phi=g(\alpha, \delta)$ basta substituir $\alpha$ e $\delta$ por seus respectivos estimadores de máxima verossimilhança. Essa é outra propriedade importante do estimador de máxima verossimilhança. Se além de estimar $\phi$ existe interesse em construir um intervalo de confiança, então é necessário obter uma estimativa para o erro padrão de $\phi.$ Isso é feito usando o Método Delta que é descrito a seguir.
Considere inicialmente que θ é um escalar e desejamos avaliar $Var[g(\widehat{\theta})].$ Expandindo $g(\widehat{\theta})$ em torno de
$$E[\widehat{\theta}] =\theta$$
e ignorando os termos superiores ao de primeira ordem temos
$$g(\widehat{\theta})=g(\theta)+(\widehat{\theta}-\theta)\dfrac{\partial g(\theta)}{\partial\theta}$$
assim, obtemos que
$$Var(g(\widehat{\theta}))=Var(\widehat{\theta})\left(\dfrac{\partial g(\theta)}{\partial(\theta)}\right)^2$$
A versão multivariada do método delta será necessária para o caso das distribuições que envolvem mais de um parâmetro. Suponha, como anteriormente, que $\theta=(\delta, \alpha)$ e que estamos interessados em $\phi=g(\alpha, \delta).$ Procedendo de forma similar temos que
$$Var(\widehat{\phi})=Var(\widehat{\alpha})\left(\dfrac{\partial \phi}{\partial \alpha}\right)^2+Cov(\widehat{\alpha},\widehat{\delta})\dfrac{\partial\phi\partial\phi}{\partial\alpha \partial\delta}+Var(\widehat{\delta})\left(\dfrac{\partial \phi}{\partial\delta}\right)^2.$$
4.3 - Escolha do Modelo Probabilístico
A escolha do modelo a ser utilizado é um tópico extremamente importante na análise paramétrica de dados de confiabilidade.
A utilização de um modelo inadequado para representar os dados compromete a análise estatística, gerando erros grosseiros de estimativas.
Figura 1.4.13: Representação de um modelo inadequado
No entanto, se um modelo probabilístico adequado aos dados é definido, podemos então utilizar o método de máxima verossimilhança para estimar os seus parâmetros.
As distribuições de probabilidade apresentadas anteriormente são típicas para dados de confiabilidade. Entretanto, para avaliar qual distribuição melhor se ajusta aos dados podemos comparar os valores estimados com os valores observados. Portanto, a forma mais simples de avaliar qual distribuição será utilizada para um conjunto de dados e escolher um modelo mais adequado é através de técnicas gráficas.
Por fim, se nenhum modelo paramétrico for adequado, necessita-se de uma análise por meio de técnicas não-paramétricas, como o estimador de Kaplan-Meier já discutido anteriormente. No entanto, se a distribuição de probabilidade for bem especificada, as técnicas paramétricas serão mais eficientes do que as técnicas não-paramétricas.
A seguir, apresentamos as técnicas gráficas utilizadas para a escolha de um modelo probabilístico.
4.3.1 - Métodos Gráficos
Método 1: Comparação da Função de Confiabilidade do Modelo Proposto com o Estimador de Kaplan-Meier
Aqui ajustamos os modelos propostos ao conjunto de dados (digamos, os modelos lognormal e de Weibull) e a partir das estimativas dos parâmetros de cada modelo estimamos a função de confiabilidade. Comparamos graficamente a função de confiabilidade de cada modelo proposto com a função de confiabilidade de Kaplan-Meier. Assim, o “melhor” modelo é aquele cujos pontos da função de confiabilidade estão mais próximos dos valores obtidos pela estimativa de Kaplan-Meier.
Método 2: Linearização da função de Confiabilidade - Papel de Probabilidade
Nesse caso, a ideia básica consiste em construir gráficos que devem ser aproximadamente lineares caso o modelo proposto seja apropriado. Violações da linearidade podem ser verificadas visualmente. A seguir, são apresentadas as técnicas de linearização para os modelos de Weibull e Log-normal.
Modelos Weibull e Exponencial
A partir da equação de confiabilidade da distribuição de Weibull, obtemos a seguinte relação
$$\log[-\log(R(t))]=\delta \log(t) - \delta \log(\alpha).$$
Portanto, uma maneira de verificar a adequação da distribuição de Weibull aos dados é plotar o gráfico de $\log[-\log(\widehat{R}(t)]$ contra log(t) e verificar se ele é aproximadamente linear, sendo $\widehat{R}(t)$ a estimativa de Kaplan-Meier para a função de confiabilidade.
Modelo Log-Normal
Analogamente, a partir da equação de confiabilidade da distribuição log-normal concluímos que
$$\Phi^{-1}[R(t)]= -\dfrac{1}{\sigma} \log(t) + \dfrac{\mu}{\sigma}, \tag{4.3.1.1}$$
sendo que $\Phi^{-1}(\cdot)$ é uma função que associa a cada 0 < p < 1 o quantil 100×p% da distribuição Normal padrão. Os valores de $\Phi^{-1}(p)$ para cada p podem ser obtidos a partir de uma tabela ou por meio de algum software estatístico.
Da equação (4.3.1.1) concluímos que se o modelo log-normal for adequado então o gráfico $\Phi^{-1}[\widehat{R}(t)]$ contra $\log(t)$ deve ser aproximadamente linear.
Observe que, a partir dos dois gráficos é possível obter estimativas grosseiras para os parâmetros dos modelos. Por exemplo, no caso do modelo de Weibull podemos encontrar uma reta que melhor se ajusta aos pontos do gráfico $\log[-\log(\widehat{R} (t))]$ contra $\log(t)$. A inclinação dessa reta é uma estimativa para δ e o intercepto uma estimativa para $\delta \log(\alpha)$. Analogamente, para o caso do modelo log-normal a inclinação do gráfico $\Phi^{-1}[\widehat{R}(t)]$ contra $\log(t)$ é uma estimativa para $-1/σ$ e o intercepto uma estimativa para $μ/σ$.
A forma mais apropriada para se obter estimativas dos parâmetros após selecionar o modelo é aplicar o método da máxima verossimilhança. Entretanto, em muitos casos a obtenção das estimativas de máxima verossimilhança dependem da utilização de métodos computacionais que funcionam melhor quando se fornecem valores iniciais para os parâmetros. Com isso, é usual a obtenção desses valores iniciais a partir de estimativas baseadas na reta de regressão de $\log[-\log(\widehat{R}(t))]$ contra $\log(t)$.
4.3.2 - Teste de Qualidade de Ajuste de Anderson Darling
Uma outra forma de checar a adequação do modelo probabilístico aos dados é testando a hipótese de que uma dada amostra tenha sido retirada de uma população com função de distribuição acumulada (f.d.a.) contínua F(x).
Seja $x_{1}$, $x_{2}$, …, $x_{n}$ uma amostra aleatória e suponha que um provável candidato para a f.d.a. dos dados seja F(x), o teste de hipóteses para verificar a adequação da distribuição $F(x)$ aos dados é
$$\begin{cases} H_{0}: \hbox{a amostra tem distribuição~} F(x), \cr H_{1}: \hbox{a amostra não tem distribuição~} F(x). \end{cases} \tag{4.3.2.1}$$
Anderson e Darling (1952, 1954) propuseram a seguinte estatística para testar (4.3.2.1)
$$A^{2} = n~\int_{-\infty}^{\infty} \dfrac{\left[ F_{n}(x)-F(x)\right] }{F(x)(1-F(x))}dF(x),$$
sendo $F_{n}(x)$ a função de distribuição acumulada empírica definida por
$$F_{n}(x)=\dfrac{1}{n}\sum_{i=1}^{n} \mathbb{I}_{(x_i\leq x)}=\begin{cases}0,~~\hbox{se}~x~<~x_{(1)}, \cr \cr \dfrac{k}{n},~~\hbox{se}~x_{(k)} \leq x~<~x_{(k + 1)}, \cr \cr 1,~~\hbox{se}~x~>~x_{(n)},\end{cases}$$
sendo $x_{(1)}\leq x_{(2)}\leq …\leq x_{(n)},$ as estatísticas de ordem da amostra aleatória e $\mathbb{I}_{(x_i \leq x)}$ a função indicadora que vale 1 se $x_i \leq x$ e 0 se $x_i~>~x$, i = 1, …, n.
A estatística $A^2$ pode ser representada numa forma equivalente como
$$A^{2}=-n-\dfrac{1}{n} \sum^{n}_{i=1}\left[(2i-1)\ln(F(x_{(i)}))+(2(n-i)+1)\ln(1-F(x_{(i)}))\right]$$
A transformação $F(x_{(i)})$ leva $x_{(i)}$ em $U_{(i)}$, sendo $U_{(1)}$, …, $U_{(n)}$ uma amostra de tamanho n com distribuição uniforme em (0,1). Logo,
$$A^{2} = - n - \dfrac{1}{n}\sum^{n}_{i=1}\left[(2i - 1)\ln(U_{(i)})+ (2(n-i) + 1)\ln(1 - U_{(i)}) \right] \qquad \tag{4.3.2.2}$$
Para calcular o valor da estatística $A^2$, devemos seguir os passos abaixo:
-
Ordene os valores da amostra: $x_{(1) }$≤ $x_{(2) }$≤ … ≤ $x_{(n)}$;
-
Quando necessário, estime os parâmetros da distribuição de interesse;
-
Calcule $U_{i}$ = $F(x_{(i)})$ e calcule o valor da estatística de Anderson Darling (4.3.2.2):
$$A^{2} = - n - \dfrac{1}{n}\sum^{n}_{i=1}\left[(2i - 1)\left(\ln(~U_{i})+ \ln(1 - U_{n+1-i}~)\right) \right]$$
(observe que esta é uma forma equivalente à (4.3.2.2))
- Para cada uma das distribuições calcule, se for o caso, o valor da estatística modificada de acordo com as tabelas dadas para cada uma delas.
Para uma distribuição com parâmetros conhecidos podemos encontrar os valores da função de distribuição acumulada da estatística $A^2$ tabulados em Peter and Lewis(1960). No entanto, surge um problema quando um ou dois dos parâmetros da distribuição precisam ser estimados, para contornar esse problema Stephens(1974, 1976, 1977) utilizou métodos assintóticos para tabular os valores dessas probabilidades quando os parâmetros das distribuições são desconhecidos.
Aplicação
Vamos aplicar o teste de qualidade de ajuste de Anderson Darling a algumas das distribuições de probabilidade mais conhecidas tais como a Normal, Exponencial, Weibull, Lognormal e Valor Extremo. Para essas distribuições o parâmetro θ pode ser univariado ou bivariado, isto é, ele terá no máximo duas componentes, conforme os seguintes casos:
-
Caso 0: O parâmetro $\theta = (\alpha, \beta)$ é totalmente conhecido;
-
Caso 1: $\alpha$ é conhecido;
-
Caso 2: $\beta$ é conhecido;
-
Caso 3: Nenhum dos componentes de $\theta= (\alpha, \beta)$ é conhecido.
Distribuição Normal
Para a distribuição Normal com função densidade de probabilidade dada por
$$f(x) = \dfrac{1}{\sqrt{2\pi}\sigma} \exp\left(-\dfrac{(x-\mu)^{2}}{2 \sigma^{2}}\right), \qquad -\infty~<~x~<~\infty.$$
A Tabela 1.4.1 fornece alguns valores para os quantis da distribuição da estatística de Anderson Darling modificada de acordo com cada um dos casos:
-
Caso 0: O parâmetro $θ = (μ, σ)$ é totalmente conhecido;
-
Caso 1: $μ$ é conhecido e $σ$ é estimado por $s^2$;
-
Caso 2: $σ$ é conhecido e $μ$ é estimado por $\overline{x}$;
-
Caso 3: Nenhum dos componentes de $θ = (μ, σ)$ é conhecido e são estimados por $(\overline{x}, s^{2})$.
A Tabela 1.4.1 fornece os quantis $1- \alpha$ da distribuição de $A^2$, ou seja, fornecem pontos $q_{(1-\alpha)}$ para os quais a probabilidade de $A^2$ ser maior que $q_{(1-\alpha)}$ é igual a $\alpha$.
| Pontos | percentis | para | cada | $\alpha$ (%) | ||
|---|---|---|---|---|---|---|
| Caso | Modificação | 15,0 | 10,0 | 5,0 | 2,5 | 1,0 |
| 0 | Nenhuma | 1,610 | 1,933 | 2,492 | 3,070 | 3,857 |
| 1 | - | 0,784 | 0,897 | 1,088 | 1,281 | 1,541 |
| 2 | - | 1,433 | 1,761 | 2,315 | 2,890 | 3,682 |
| 3 | $A^2(1+(4/n)-(25/n^2))$ | 0,560 | 0,632 | 0,751 | 0,870 | 1,029 |
Tabela 1.4.1: Tabela de pontos percentis de $A^2$ para a distribuição Normal.
Exemplo 4.3.2.1
Considere as seguintes medidas correspondentes ao peso de homens (em pounds): 148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236.
Desejamos testar a seguinte hipótese:
$$\begin{cases}H_{0}: \hbox{Os dados seguem distribuição Normal~} - N(\mu, \sigma) \cr H_{1}: \hbox{Os dados não seguem distribuição Normal~} - N(\mu, \sigma)\end{cases}$$
A média dos dados é $\overline{x} = 172$ e o desvio padrão é $s = 24,9520$.
| dados | dados ordenados | $F_{(xi)}$ | $ln(F_{(xi)})$ | $ln(1 - F_{(xi)})$ |
|---|---|---|---|---|
| 148 | 148 | 0,168063 | -1,78341 | -0,184 |
| 154 | 154 | 0,235336 | -1,44674 | -0,26832 |
| 158 | 158 | 0,287372 | 1,24698 | -0,3388 |
| 160 | 160 | 0,315285 | -1,15428 | 0,37875 |
| 161 | 161 | 0,329662 | -1,10969 | -0,39997 |
| 162 | 162 | 0,344295 | -1,06626 | -0,42204 |
| 166 | 166 | 0,404986 | -0,9039 | -0,51917 |
| 170 | 170 | 0,468057 | -0,75916 | -0,63122 |
| 182 | 182 | 0,655705 | -0,42204 | -1,06626 |
| 195 | 195 | 0,821676 | -0,19641 | -1,72415 |
| 236 | 236 | 0,99484 | -0,00517 | -5,26684 |
Tabela 1.4.2: Calculando o valor de $A^2$.
Utilizando a fórmula (4.3.2.2), temos que
$$ D = (2×1 - 1)×(-1,78341) + (2×(11 - 1) + 1)×(-0,184)$$ $$ + (2×2 - 1)×(-1,44674) + (2×(11 - 2) + 1)×(-0,26832)$$ $$+ (2×3 - 1)×(-1,24698) + (2×(11 - 3) + 1)×(-0,3388)$$ $$+ (2×4 - 1)×(-1,15428) + (2×(11 - 4) + 1)×(-0,37875) $$ $$+ (2×5 - 1)×(-1,10969) + (2×(11 - 5) + 1)×(-0,39997) $$ $$+ (2×6 - 1)×(-1,06626) + (2×(11 - 6) + 1)×(-0,42204)$$ $$+ (2×7 - 1)×(-0,9039) + (2×(11 - 7) + 1)×(-0,51917)$$ $$+ (2×8 - 1)×(-0,75916) + (2×(11 - 8) + 1)×(-0,63122)$$ $$+ (2×9 - 1)×(-0,42204) + (2×(11 - 9) + 1)×(-1,06626)$$ $$+ (2×10 - 1)×(-0,19641) + (2×(11-10) + 1)×(-1,72415)$$ $$+ (2×11 - 1)×(-0,00517) + (2×(11-11) + 1)×(-5,26684)$$ $$ = -131,4145$$
Com isso, temos
$$A^{2} = -\dfrac{D}{n} - n = \dfrac{131,4145}{11} - 11 = 0,9467719.$$
A estatística de Anderson Darling modificada para o Caso 3 (μ e σ desconhecidos) é dada por:
$$A^{2}_{m} = A^{2} \times (1 + (4/n) - (25/n^{2}))= A^{2} \times (1 + (4/n) - (25/n^{2}))=0,94677 \times 1,15703 = 1,0954.$$
A partir da Tabela 1.4.1 concluímos que o p-valor do teste é menor que 0,01. Portanto, assumindo um nível de significância igual a 0,05 rejeitamos a hipótese dos dados serem provenientes de uma distribuição normal.
Distribuição Log-Normal
Para realizar o teste de Anderson Darling quando a distribuição dada é Log-normal devemos considerar o logaritmo dos dados e proceder como no caso da distribuição Normal.
Distribuição Exponencial
Considere o teste como em (4.3.2.1) e a distribuição Exponencial com função de distribuição acumulada dada por
$$F(x) = 1- \exp\left(-\dfrac{x}{\alpha}\right), \qquad x~>~0.$$
Os seguintes casos podem ocorrer durante a realização do teste:
-
Caso 0: O parâmetros $\alpha$ é conhecido;
-
Caso 1: O parâmetro $\alpha$ precisa ser estimado.
A Tabela 1.4.3 apresenta os valores da estatística $A^2$ com modificações apropriadas para cada um dos casos citados acima.
| Pontos | percentis | para | cada | $\alpha$ (%) | ||
|---|---|---|---|---|---|---|
| Caso | Modificação | 15,0 | 10,0 | 5,0 | 2,5 | 1,0 |
| 0 | Nenhuma | 1,610 | 1,933 | 2,492 | 3,070 | 3,857 |
| 3 | $A^2(1+(0,6/n))$ | 0,922 | 1,078 | 1,341 | 1,606 | 1,957 |
Tabela 1.4.3: Tabela de pontos percentis de $A^2$ para a distribuição Exponencial.
Exemplo 4.3.2.2
Os dados a seguir se referem aos tempos de vida (em horas) de 15 componentes eletrônicos colocados em teste. Sejam eles: 7,134; 1,157; 103,507; 64,707; 48,826; 72,332; 155,894; 83,653; 5,729; 4,472; 14,578; 42,833; 45,118; 223,395; 3,055. A média dos dados é $\overline{x}=58.4260$.
Desejamos realizar o seguinte teste:
$$\begin{cases}H_{0}:~\hbox{Os dados seguem distribuição Exponencial~}-\hbox{Exp}(\lambda) \cr H_{1}: \hbox{Os dados não seguem uma distribuição Exponencial~}- \hbox{Exp}(\lambda) \end{cases}$$
Assim, para calcular o valor da estatística $A^2$ procedemos como na Tabela 1.4.4.
| dados | dados ordenados | $F_{(xi)}$ | $ln(F_{(xi)})$ | $ln(1 - F_{(xi)})$ |
|---|---|---|---|---|
| 7,134 | 1,157 | 0,019603 | -3,93208 | -0,0198 |
| 1,157 | 3,055 | 0,050948 | -2,97696 | -0,05229 |
| 103,507 | 4,472 | 0,073685 | -2,60796 | -0,07654 |
| 64,707 | 5,729 | 0,0934 | -2,37086 | -0,09805 |
| 48,826 | 7,134 | 0,114939 | -2,16336 | -0,1221 |
| 72,332 | 14,578 | 0,220821 | -1,5104 | -0,24952 |
| 155,894 | 42,833 | 0,519593 | -0,65471 | -0,73312 |
| 83,653 | 45,118 | 0,538019 | -0,61986 | -0,77223 |
| 5,729 | 48,826 | 0,566423 | -0,56841 | -0,83569 |
| 4,472 | 64,707 | 0,669617 | -0,40105 | -1,1075 |
| 14,578 | 72,332 | 0,710041 | -0,34243 | -1,23802 |
| 42,833 | 83,653 | 0,761115 | -0,27297 | -1,43177 |
| 45,118 | 103,507 | 0,829938 | -0,1864 | -1,77159 |
| 223,395 | 155,894 | 0,930625 | -0,0719 | -2,66823 |
| 3,055 | 223,395 | 0,97815 | -0,02209 | -3,82356 |
Tabela 1.4.4: Calculando o valor de $A^2$.
Com isso, temos
$$D = (2×1 - 1)×(-3,93208) + (2×(15 - 1) + 1)(-0,01980)$$ $$ + (2×2 - 1)×(-2,97696) + (2×(15 - 2) + 1)(-0,05229)$$ $$ + (2×3 - 1)×(-2,60796) + (2×(15 - 3) + 1)(-0,07654)$$ $$ + (2×4 - 1)×(-2,37086) + (2×(15 - 4) + 1)(-0,09805)$$ $$ + (2×5 - 1)×(-2,16336) + (2×(15 - 5) + 1)(-0,12210)$$ $$ + (2×6 - 1)×(-1,51040) + (2×(15 - 6) + 1)(-0,24952)$$ $$ + (2×7 - 1)×(-0,65471) + (2×(15 - 7) + 1)(-0,73312)$$ $$ + (2×8 - 1)×(-0,61986) + (2×(15 - 8) + 1)(-0,77223)$$ $$ + (2×9 - 1)×(-0,56841) + (2×(15 - 9) + 1)(-0,83569)$$ $$ + (2×10 - 1)×(-0,40105) + (2×(15 - 10) + 1)(-1,10750)$$ $$ + (2×11 - 1)×(-0,34243) + (2×(15 - 11) + 1)(-1,23802)$$ $$ + (2×12 - 1)×(-0,27297) + (2×(15 - 12) + 1)(-1,43177)$$ $$ + (2×13 - 1)×(-0,18640) + (2×(15 - 13) + 1)(-1,77159)$$ $$ + (2×14 - 1)×(-0,07190) + (2×(15 - 14) + 1)(-2,66823)$$ $$ + (2×15 - 1)×(-0,02209) + (2×(15 - 15) + 1)(-3,82356)$$ $$ = -236,79011 $$
Assim, temos que
$$A^{2} = -\dfrac{D}{n} - n = -\dfrac{-236,79011}{15} - 15 = 0,786007333$$
Logo, a estatística de Anderson Darling modificada de acordo com a Tabela 1.4.3 é dada por:
$$A^{2}_{m} = A^{2}(1 + (0,6/n)) = 0,786007333\times(1+(0,6/15)) = 0,8174$$
Com isso, a partir da Tabela 1.4.3 concluímos que o p-valor é maior que 0,15.
Distribuição Valor Extremo
Para realizar o teste (4.3.2.1) para a distribuição de Valor Extremo com função de distribuição acumulada dada por
$$F(y) = 1 - \exp \left[-\exp\left(\dfrac{y-\mu}{\sigma}\right)\right], \qquad -\infty~<~y~<~\infty, \tag{4.3.2.3}$$
Os seguintes casos podem ocorrer:
-
Caso 0: O parâmetro θ = (μ, σ) é totalmente conhecido;
-
Caso 1: O parâmetro μ é conhecido e σ precisa ser estimado;
-
Caso 2: O parâmetro σ é conhecido e μ precisa ser estimado;
-
Caso 3: Nenhum dos componentes de θ = (μ, σ) é conhecido e portanto ambos precisam ser estimados.
A seguinte tabela apresenta os valores da estatística $A^2$ com modificações apropriadas para cada um dos casos citados acima.
| Pontos | percentis | para | cada | $\alpha$ (%) | ||
|---|---|---|---|---|---|---|
| Caso | Modificação | 15,0 | 10,0 | 5,0 | 2,5 | 1,0 |
| 0 | Nenhuma | - | 1,933 | 2,492 | 3,070 | 3,857 |
| 1 | $A^2(1+(0,3/n))$ | 0,736 | 1,062 | 1,321 | 1,591 | 1,959 |
| 2 | Nenhuma | 1,060 | 1,725 | 2,277 | 2,854 | 3,640 |
| 3 | $A^2(1+(0,2/\sqrt{n}))$ | 0,474 | 0,637 | 0,757 | 0,877 | 1,038 |
Tabela 1.4.5: Tabela de pontos percentis de $A^2$ para a Valor Extremo.
Exemplo 4.3.2.3
Considere os dados a seguir provenientes de uma distribuição de Valor Extremo: 84,01; 75,498; 79,356; 72,635; 104,052; 102,56; 91,458; 90,546; 78,932; 90,18; 76,828; 93,905; 75,433; 85,35; 102,64.
Os parâmetros estimados são: locação $μ = 92,162789$ e escala $σ = 9,976448$.
O objetivo é realizar o seguinte teste:
$$\begin{cases}H_{0}: \hbox{Os dados seguem distribuição Valor Extremo~}- \hbox{VE}(\mu,\sigma) \cr H_{1}:~\hbox{Os dados não seguem distribuição Valor Extremo~}- \hbox{VE}(\mu,\sigma) \end{cases}$$
Assim, vamos calcular o valor da estatística $A^2$ procedendo como na Tabela 1.4.6.
| dados | dados ordenados | $F_{(xi)}$ | $ln(F_{(xi)})$ | $ln(1 - F_{(xi)})$ |
|---|---|---|---|---|
| 84,01 | 72,635 | 0,131708 | -2,02717 | -0,14123 |
| 75,498 | 75,433 | 0,170513 | -1,76895 | -0,18695 |
| 79,356 | 75,498 | 0,171526 | -1,76302 | -0,18817 |
| 72,635 | 76,828 | 0,193462 | -1,64268 | -0,21500 |
| 104,052 | 78,932 | 0,233165 | -1,45601 | -0,26548 |
| 102,56 | 79,356 | 0,241953 | -1,41901 | -0,27701 |
| 91,458 | 84,01 | 0,357035 | -1,02992 | -0,44167 |
| 90,546 | 85,35 | 0,396589 | -0,92486 | -0,50516 |
| 78,932 | 90,18 | 0,559461 | -0,58078 | -0,81976 |
| 90,18 | 90,546 | 0,572752 | -0,55730 | -0,85039 |
| 76,828 | 91,458 | 0,606153 | -0,50062 | -0,93179 |
| 93,905 | 93,905 | 0,696025 | -0,36237 | -1,19081 |
| 75,433 | 102,56 | 0,941304 | -0,06049 | -2,83538 |
| 85,35 | 102,64 | 0,942629 | -0,05908 | -2,85821 |
| 102,64 | 104,052 | 0,962849 | -0,03786 | -3,29277 |
Tabela 1.4.6: Calculando o valor de $A^2$.
Assim, temos que
$$A^{2} = 0,6421041$$
Logo, a estatística de Anderson Darling modificada de acordo com a Tabela 1.4.5 é dada por:
$$A^{2}_{m} = A^{2}(1 + (0,2/\sqrt{15}))= 0,6421041\times 1,051640= 0,6753.$$
A partir da Tabela 1.4.5, concluímos que o p-valor está entre 0,05 e 0,10. Para obtermos um valor exato, podemos fazer uma interpolação entre esses valores a partir da equação que fornece a inclinação da reta, dada por
$$\dfrac{0,757 - 0,637}{5 - 10} = \dfrac{0,6753 - 0,637}{x - 10}.$$
Dessa equação, concluímos que o p-valor é de 8,4042% ou 0,0840. Portanto, para um nível de significância igual a 0,05, não rejeitamos a hipótese $H_{0}$.
Distribuição Weibull
Para realizar o teste (4.3.2.1) para a distribuição de Weibull, tomamos o logaritmo dos dados e procedemos como no caso da distribuição de Valor Extremo. Por exemplo, se a variável aleatória X tem distribuição Weibull com função distribuição acumulada dada por:
$$F(x) = 1 - \exp\left(-\left(\dfrac{x}{\alpha}\right)^{\delta}\right), \qquad x~>~0,$$
então a variável aleatória $Y = \log(X)$ tem distribuição Valor Extremo, dada por (4.3.2.3), com parâmetros de locação $μ$ e de escala $σ$ ($\sigma = 1/\delta$ e $\mu = \log(\alpha)$).
4.3.3 - Aplicação
Exemplo 4.3.3.1
O fabricante de um tipo de isolador elétrico quer conhecer o comportamento de seu produto funcionando na temperatura de 150ºC. Um teste de vida foi realizado nestas condições usando 45 isoladores elétricos. O teste terminou quando 30 deles haviam falhado (censura do tipo II). As 15 unidades restantes que não haviam falhado foram censuradas no tempo t = 836. Os tempos das falhas (em horas) estão dispostos na tabela abaixo, onde o símbolo + denota censura.
| 123 | 230 | 250 | 261 | 276 | 291 | 296 | 299 | 311 | 332 |
| 338 | 379 | 387 | 425 | 448 | 472 | 479 | 488 | 489 | 502 |
| 565 | 589 | 603 | 603 | 605 | 626 | 635 | 660 | 661 | 836 |
| 836+ | 836+ | 836+ | 836+ | 836+ | 836+ | 836+ | 836+ | 836+ | 836+ |
| 836+ | 836+ | 836+ | 836+ | 836+ |
Tabela 1.4.7: Tempos das falhas (em horas)
O fabricante tem interesse em estimar o tempo médio (MTTF) de vida do isolador e o percentual de falhas até 300 horas de uso.
Para isso, devemos seguir os passos abaixo:
1. Escolher o modelo probabilístico: vamos utilizar o Método 2, estudado no tópico anterior.
O Gráfico QQ-plot para as distribuições Normal, Exponencial, Weibull e Log-normal é apresentado na Figura 1.4.14.
Figura 1.4.14: Gráficos QQ-plot.
O valor da estatística e do p-valor do teste de Anderson-Darling é dado na Tabela 1.4.8. Podemos observar que para a distribuição Weibull, o valor da estatística de Anderson-Darling é o menor dentre todos. Logo, concluímos que a distribuição de Weibull é a mais adequada para esse conjunto de dados.
| Distribuição | Anderson-Darling | P-valor |
|---|---|---|
| Normal(mu = 448.63, sigma = 170) | 0.4205505153 | 0.30453 |
| Exponencial(Taxa = 0.00222899) | 5.4881391164 | 0.00001 |
| Weibull(Forma = 3.01512, Escala = 503.113) | 0.3848473211 | 0.25 |
| Log-Normal Log-Normal(log(mu) = 6.03012, log(sigma) = 0.41011) | 0.5312456754 | 0.16054 |
Tabela 1.4.8: Valor da estatística e do p-valor do teste de Anderson-Darling.
2. Análise do tempo de falha utilizando a distribuição de Weibull
Usando o software Action e o método de estimação de máxima verossimilhança obtemos as seguintes estimativas dos parâmetros:
| Parâmetros | Estimativas | Desvio-padrão | Lim. Inferior | Lim. Superior |
|---|---|---|---|---|
| Forma (δ) | 1,9361 | 0,0943 | 1,7512 | 2,1209 |
| Escala (α) | 762,0024 | 5283,0479 | 0,0000 | 11116,5860 |
Tabela 1.4.9: Tabela de estimativas da distribuição Weibull
Com os valores da tabela acima, a estimativa de máxima verossimilhança para o MTTF é dada por
$$\widehat{MTTF}=\widehat{\alpha} , \Gamma\left[1 + \dfrac{1}{\widehat{\delta}}\right]$$
$$=762,0024\times\Gamma\left[1+\dfrac{1}{1,9361} \right]=762,0024\times\Gamma\left[1,5165 \right]$$
$$=762,0024 \times 0,8869 = 675,8199.$$
Analogamente, a estimativa de máxima verossimilhança para a função de confiabilidade R(t) é dada por
$$\widehat{R}(t)=\exp\left[-\left(\dfrac{t}{\widehat{\alpha}}\right)^{\widehat{\delta}}\right]=\exp\left[-\left(\dfrac{t}{762,0024}\right)^{1,9361}\right].$$
Portanto,
$$\widehat{R}(300)=\exp\left[-\left(\dfrac{300}{762,0024}\right)^{1,9361}\right]$$
$$=\exp\left[-\left(0,3937\right)^{1,9361}\right]=\exp\left[-0,1645\right]=0,8483,$$
isto é, estima-se que a probabilidade de um isolador elétrico operar por mais de 300 horas é de 84,83%. Portanto, a probabilidade de um isolador falhar antes de 300 horas é de aproximadamente 15%.
Exercício 4.3.3.1
No Exemplo 2.1, envolvendo um teste de vida com as válvulas, o fabricante estava interessado em estimar as seguintes características de seu produto:
-
o tempo médio de falha (MTTF);
-
o tempo no qual 10% dos mecanismos estão fora de operação;
-
o percentual esperado de falha nos dois primeiros anos de uso.
Escolha um modelo probabilístico adequado e responda às questões do fornecedor.
4.4 - Método Delta
Muitas vezes temos interesse em calcular intervalos de confiança para funções dos parâmetros de uma determinada distribuição. Por exemplo, no caso de uma distribuição Weibull, intervalos de confiança para o MTTF, intervalos de confiança para os percentis ou para a mediana. Todos essas quantidades são do tipo $\phi=g(\alpha,\delta)$ e para estimar $\phi$ basta tomar os estimadores de máxima verossimilhança e substituir em $\phi$, obtendo $\widehat{\phi}=g(\widehat{\alpha},\widehat{\delta})$. Essa é a propriedade de invariância conhecida dos estimadores de máxima verossimilhança. Para calcular o intervalo de confiança para $\phi=g(\alpha,\delta)$ teremos que calcular o erro padrão de $\widehat{\phi}=g(\widehat{\alpha},\widehat{\delta})$.
Tomando $g(\widehat{\alpha},\widehat{\delta})$ e expandindo em série de Taylor até a primeira ordem, temos
$$g(\widehat{\alpha},\widehat{\delta}) \approx g(\alpha,\delta)+(\widehat{\alpha}-\alpha) \left(\dfrac{\partial g(\alpha,\delta)}{\partial \alpha}\right)+(\widehat{\delta}-\delta) \left(\dfrac{\partial g(\alpha,\delta)}{\partial\delta}\right)+2(\widehat{\alpha}-\alpha)(\widehat{\delta}-\delta)\left(\dfrac{\partial^2 g(\alpha,\delta)}{\partial \alpha~\partial\delta}\right).$$
Assim, temos
$$Var\left(g(\widehat{\alpha},\widehat{\delta})\right)\approx Var\left(\widehat{\alpha}\right)\left(\dfrac{\partial g(\alpha,\widehat{\delta})}{\partial \alpha}\right)^2 + Var\left(\widehat{\delta}\right)\left(\dfrac{\partial g(\alpha,\delta)}{\partial \delta}\right)^2 + 2~Cov\left(\widehat{\alpha},\widehat{\delta}\right)\left(\dfrac{\partial^2g(\alpha,\delta)}{\partial \alpha ~\partial \delta}\right)$$
$$\approx Var\left(\widehat{\alpha}\right)\left(\dfrac{\partial \widehat{\phi}}{\partial \alpha}\right)^2 + Var\left(\widehat{\delta}\right)\left(\dfrac{\partial \widehat{\phi}}{\partial \delta}\right)^2 + 2~Cov\left(\widehat{\alpha},\widehat{\delta}\right)\left(\dfrac{\partial \widehat{\phi}}{\partial \alpha}\dfrac{\partial \widehat{\phi}}{\partial \delta}\right)$$
Vamos utilizar o método delta para calcular um intervalo de confiança com nível de significância de 5% (nível de confiança 95%) para o quantil da distribuição Weibull, utilizando os dados do Exemplo 4.3.3.1. Como estimador do quantil temos a seguinte função dos parâmetros
$$\widehat{t_p} = \widehat{\alpha} \left(-\log(1 - p)\right)^{1/\widehat{\delta}}$$
Com isso,
$$\widehat{Var\left(\widehat{t_p}\right) } = Var(\widehat{\alpha}) \ast A^2 + Var(\widehat{\delta}) \ast B^2 + 2 \ast Cov(\widehat{\alpha},\widehat{\delta}) \ast A \ast B$$
em que
$$A= \left(-\log(1 - p)\right)^{1 / \widehat{\delta}}$$
$$B= \dfrac{\displaystyle \widehat{\alpha}\left(\left(-\log(1 - p)\right)^{1/\widehat{\delta}}\log\left( \left(-\log(1 - p)\right )^{1/\widehat{\delta}}\right)\right)}{\displaystyle \widehat{\delta}^2}=\dfrac{\displaystyle \widehat{\alpha} \ast A \ast \log(A)}{\displaystyle \widehat{\delta}^2}$$
4.5 - Resumo das Principais Distribuições e Propriedades
A seguir apresentamos um resumo das principais distribuições e suas propriedades.
1. Distribuição Exponencial
1.1 Função densidade de probabilidade
$$f(t) = \alpha \exp\left(-\alpha t\right),~~~~~~t~>~0$$
em que $\alpha~>~0$ é o parâmetro de escala.
1.2 Função de confiabilidade
$$R(t)=\exp\left(-\alpha t\right).$$
1.3 Função taxa de falha
$$h(t)=\alpha.$$
1.4 MTTF e variância
$$\textrm{MTTF}=\dfrac{1}{\alpha}$$
$$Var(T)=\dfrac{1}{\alpha^2}.$$
1.5 Quantil 100×p%:
$$t_p=-\dfrac{1}{\alpha} \log(1 - p).$$
1.6 Relação com a distribuição de Valor Extremo
Se T tem distribuição exponencial com parâmetro $\alpha$ então log(T) tem distribuição de Valor Extremo com parâmetros $\mu= \log(\alpha)~$ e $~\sigma=1.$
1.7 Linearização de R(t)
$$\log[-\log(R(t))]=\log(t)+\log(\alpha).$$
2. Distribuição de Weibull
2.1 Função densidade de probabilidade
$$f(t)=\dfrac{\delta}{\alpha^\delta}~t^{\delta-1}\exp\left[-\left(\dfrac{t}{\alpha}\right)^\delta\right],~~~~~~t~>~0$$
em que $\delta~>~0$ e $\alpha~>~0$ são, respectivamente, os parâmetro de forma e escala.
2.2 Função de confiabilidade
$$R(t)=\exp\left[-\left(\dfrac{t}{\alpha}\right)^\delta\right].$$
2.3 Função taxa de falha
$$h(t)=\dfrac{\delta}{\alpha^\delta}~t^{\delta-1}.$$
2.4 MTTF e variância
$$\textrm{MTTF}=\alpha , \Gamma\left[1 + \dfrac{1}{\delta}\right]$$
$$Var(T)=\alpha^2,\left(\Gamma\left[1+\dfrac{2}{\delta}\right]-\Gamma^2\left[1+\dfrac{1}{\delta}\right]\right).$$
2.5 Quantil 100×p%
$$t_p=\alpha[-\log(1- p)]^{1/\delta}.$$
2.6 Relação com a distribuição de Valor Extremo
Se T tem distribuição Weibull ($\alpha, \delta$) então log(T) tem distribuição de Valor Extremo com parâmetros $\mu = \log(\alpha )~$ e $~\sigma = \dfrac{1}{\delta}$.
2.7 Relação com a distribuição Exponencial
Se T tem distribuição de Weibull com parâmetro de escala $\alpha$ e de forma $\delta = 1$ então T tem distribuição Exponencial com parâmetro de escala $\alpha$.
2.8 Linearização de R(t)
$$\log[-\log(R(t)]=\delta \log(t)-\delta \log(\alpha).$$
3. Distribuição Log-normal
3.1 Função densidade de probabilidade
$$f(t)=\dfrac{1}{\sqrt{2 \pi}, , t \sigma} , \exp\left[-\dfrac{1}{2}\left(\dfrac{\log(t)-\mu}{\sigma} \right)^2\right],~~~~~~t~>~0.$$
em que $-\infty~<~\mu~<~\infty~$ e $~\sigma~>0.$
3.2 Função de confiabilidade
$$R(t)=\Phi \left(-\dfrac{\log(t)-\mu}{\sigma}\right)$$
3.3 Função taxa de falha
$$h(t)=\dfrac{f(t)}{R(t)},~~~~~~t~>~0.$$
3.4 MTTF e variância
$$\textrm{MTTF}=\exp\left(\mu + \dfrac{\sigma^2}{2}\right)$$
$$Var(T)=\exp(2\mu + \sigma^2) [\exp(\sigma^2)-1].$$
3.5 Quantil 100×p%
$$t_p=\exp(z_p~\sigma+\mu).$$
3.6 Relação com a distribuição Normal
Se T tem distribuição Log-normal com parâmetros $\mu$ e $\sigma$ então log(T) tem distribuição Normal com média $\mu$ e desvio-padrão $\sigma.$
3.7 Linearização de R(t)
$$\Phi^{-1}(R(t))=-\dfrac{\log(t)}{\sigma}+\dfrac{\mu}{\sigma}.$$