Análise de Sobrevivência€¦ · Análise de Sobrevivência é o estudo de indivíduos (itens observados) onde um evento bem definido (falha) ocorre depois de algum tempo (tempo

Análise de Sobrevivência

Rodrigues Z. Fazenda

Rodrigues Fazenda 1

1. CONCEITOS ESTATÍSTICOS EM SOBREVIVÊNCIA

1.1 Introdução

Análise de Sobrevivência é o estudo de indivíduos (itens observados) onde um evento bem definido (falha) ocorre depois de algum tempo (tempo de falha).

Exemplos:(i) tempo de falha de equipamentos industriais

(engenharia)(ii) tempo de sobrevida de um paciente (medicina)(iii) tempo de duração do período de desemprego ou (iii) tempo de duração do período de desemprego ou

greve (economia)Definindo algumas características:(i) as variáveis de resposta são não-negativas(ii) principalmente univariados e contínuos(iii) presença comum de censura Primeiramente, é importante ter uma definição precisa

de tempo de falha. Isto requer especificações sobre:• origem do tempo de falha• unidade de medida do tempo (calendários, tempo de

operação, milhagem, número de ciclos)• falha (mais fácil na medicina: morte, não tão fácil na

engenharia).Rodrigues Fazenda 2

Exemplo: tempo de falha de um carroQuestões a serem feitas: Quando começar a contar o

tempo? Como medir o tempo de falha? O que é falha?Nós devemos estudar melhor todas essas

especificações.

1.2 Resultados emSobrevivência

1) Descritiva vs Inferência EstatísticaEm algumas aplicações, características descritivas

simples como média simples, função de sobrevivência egráficos de probabilidades são suficientes. Em outrasaplicações, intervalos de confiança ou taxa de influênciadedeterminadasvariáveissãoexigidas.dedeterminadasvariáveissãoexigidas.

2) CensuraUm sistema pode falhar mesmo antes que todos os

itens tenham falhado num determinado tempo. Este fatotem determinadas razões.

Os itens são normalmente censurados à direita. Maspodem ainda ser censurados à esquerda ou podemosdeterminar um intervalo de censura(mais difícil deanalisar)

3) Paramétrico ou Não-ParamétricoAmbos serão vistos no disciplina. E ainda os Semi-

Paramétricos serão apresentados.Rodrigues Fazenda 3

4) Amostras Simples vs Modelos de RegressãoSe os itens pertencem à mesma população (são

similares), então uma análise de amostras simples deveser utilizada. Se os itens não são da mesma população e sesuas diferenças podem ser contadas (máquinas submetidasa pressões distintas), então estas diferenças devem serconsideradas na Análise.Variáveis para medir tais diferenças (pressão) são

denominadas variáveis explanatórias ou covariáveis.Ambos serão vistos no decorrer da disciplina.

5) Clássica vs BayesianaAmbas serão rapidamente revisadas e vistas neste

disciplina.disciplina.Modelos Bayesianos tiveram no passado a

desvantagem de que sua Análise através de Modelos deRegressão requeriam muito esforço computacional . Isto,com o tempo, foi se tornando cada vez menos importante.

Métodos bayesianos são importantes em estudos desobrevivência porque frequentemente temos informaçõesde experiências anteriores que podem usualmente seremcombinadas com os dados e incorporadas à análise.

Rodrigues Fazenda 4

1.3 Sistema Reparáveis e Não-Reparáveis

Outra importante definição, embora somente depoisserá vista neste disciplina.

Considere o experimento com um sistema reparável eos seguintes tempos de falha cumulativos:203, 286, 481, 873, 1177, 1438, 1852, 2091, 2295, 2632.Vejamos alguns itens interessantes sobre se a taxa de

falha aumenta ou diminui com o tempo.Se o sistema é tomado como não reparável, então o

tempo entre falhas é considerado.Um exemplo de análise simples na Figura 1.1 indica o

aumento da taxa de falha.Conforme o tempo passa, é mais provável que uma

máquina deste tipo venha a falhar.Para sistemas reparáveis o exemplo da análise simples

na Figura 1.2 indica taxa de falha constante em relação aotempo.

Se os tempos entre falhas são ordenados e o sistema éreparável, a Figura 1.3 indica que a taxa de falha édecrescente.

Sistemas reparáveis com processo de tempo de falhaideais não serão vistos neste disciplina.

Rodrigues Fazenda 5

t1

H(t

1)

0 100 200 300 400

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Figura 1.1 – Função de Risco Acumulada dos Tempos entre Falhas.

t2

H(t

2)

0 500 1000 1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

Figura 1.2 – Proporção de Falhas vs Tempo de Falha.

Rodrigues Fazenda 6

t3

H(t

3)

0 500 1000 1500 2000 2500

0.0

0.2

0.4

0.6

0.8

1.0

Figura 1.3 – Proporção de Falhas vs Tempo de Falha (Ordenados)

Rodrigues Fazenda 7

1.4 Componentes e Sistemas Reparáveis

Neste disciplina, trataremos de componentesreparáveis sem referência aos sistemas que podem contertais componentes. Porém, é importante afirmar queSistemas Reparáveis também são uma importante área deestudo.

Os sistemas mais simples são:• Sistemas em Série: o sistema só funciona se todas as

componentes funcionarem.• Sistemas Paralelos: o sistema só falha se todas as

componentes falharem• Sistemask out of n : o sistemasó funcionase pelo• Sistemask out of n : o sistemasó funcionase pelo

menos k das n componentes funcionarem(se k=1→paralelo e se k=n→ série).

Esses sistemas formam grande parte dos grupos de sistemas chamados Sistemas Coerentes.

Para entender o conceito de sistemas coerentes é útildefinir o indicador de funcionamento xi para acomponente i e a função de estrutura das n componentesdo sistema:

1, se o sistema funciona.φ(x1,x2,...,xn) =

0, caso contrário.Rodrigues Fazenda 8

Sistemas coerentes têm função de estrutura quesatisfaz:

(i) φ(1,1,..., 1)=1(ii) φ(0, 0,..., 0)=0(iii) φ é não-decrescente nesses argumentos.Outros sistemas são:• Sistemas multi-estados: onde as componentes podem

estar em vários estados (não só funcionando ou falhando)• Sistemas load-sharing: onde a carga do sistema é

distribuída entre as componentes que funcionam.

1.5 Distribuições Binomial e Hipergeométrica

Análises estatísticas simples são obtidas se os temposde falha são dicotomizados:

• funcionar até certo tempo (digamos tempo de falha):defeituoso.

• funcionar além deste tempo: não-defeituoso.Distribuições Binomial e Geométrica são para um

númeroX de itens defeituosos numa amostra de tamanhon e probabilidadep do item ser defeituoso.

Se uma amostra com reposição (ou que não seja deuma população muito grande), a distribuição Binomial éobtida como:

Rodrigues Fazenda 9

P(X=k) = n pk(1-p)n-k onde 0≤ k≤ nk

E(X)=np e Var(X)=np(1-p)

Quando p é desconhecido, podemos estimá-lo como:

P = X/n e Var(X) = np(1-p)

Para n grande, np>5 e n(1-p)>5, a Binomial éaproximada pela distribuição Normal com momentosconforme os descritos anteriormente.

Para uma aproximação ao nível de significância 100(1-α)%, o intervalo de confiança para p é dado por:

(X/n - zα/2 (X(n-X)/n3)1/2, X/n + zα/2 (X(n-X)/n3)1/2)

onde zα/2 é o quartil(1-α/2) da distribuição N(0,1).Outra aproximação para n grande eλ=np (p pequeno)

é a distribuição de Poisson com médiaλ.Equivalentemente testando a independência e

constância de p.Se as amostras são sem reposição e de uma população

finita, a distribuição Hipergeométrica é obtida como:n N-n

P(X=k)= k K-k k=0,1,..,KNK

Rodrigues Fazenda 10

onde N é o tamanho da população e K é a população deitens defeituosos.

E(X)=np e Var(X)=np(1-p)N-n para p=K/NN-1

1.6 Processos de Poisson

Usado particularmente para sistemas reparáveis.Assume-se primeiramente que é observada uma série

de ocorrências em linha. (As ocorrências devem ser falhassucessivas do sistema e a linha representa o tempo real).

Assume-se que:(i) as falhas ocorrem em intervalos disjuntos e(i) as falhas ocorrem em intervalos disjuntos e

independentes(ii) ocorrência = falha(iii) a taxa de falha é uma constanteλ.

Então se X é o número de falhas num intervalo detamanho s, X tem distribuição de Poisson com médiaλs.

Também, os tempos entre falhas são independentes eexponencialmente distribuídos com MTFBλ-1.

Isto indica que a exponencial com linha base (solo) é adistribuição para o tempo de falha.

Isto pode ser generalizado para permitir taxas de falhanão constantes.

→ Processo de Poisson não-homogêneo.Rodrigues Fazenda 11

2. DISTRIBUIÇÕES DE PROBABILIDADE

2.1 Introdução

Em muitas áreas de aplicação da estatística, o pontoinicial para avaliação da variável de interesse é adistribuição Normal. Isto pode resultar de umaconsideração pragmática pura ou da argumentaçãobaseada no Teoria do Limite central, que diz que se umavariável aleatória é a soma de um grande número depequenos efeitos , então a distribuição é aproximadamenteNormal.

No contexto de sobrevivência, o caso da normalidade émuito menos usado. Para que possamos entender, temposde vida e resistênciasãoquantidadespositivas. Do pontode vida e resistênciasãoquantidadespositivas. Do pontode vista do modelo, é natural começar a pensar noprocesso de Poisson, idéias já discutidas na sessão 1.6,baseado na distribuição Exponencial. Contudo estadistribuição tem uma limitada aplicabilidade na prática,generalizações da Exponencial como a Gamma e aWeibull já provarão ter maior valor prático em modelosde sobrevivência. Estas e outras distribuições deprobabilidade comumente encontradas em estudos desobrevivência são discutidas nas sessões 2.3 à 2.7.

Outros aspectos centrais da discussão sobre análise desobrevivência são as funções de sobrevivência e de risco,e a natural ocorrência de dados censurados. Estes assuntossão discutidos nas sessões 2.2 e 2.8.Rodrigues Fazenda 12

Finalmente neste capítulo colocamos os resultadosprobabilísticos em contexto de análise de dados. Contudométodos gerais para ajustar distribuições deprobabilidades são desenvolvidos no Capítulo 3, algumastécnicas básicas são apresentadas na sessão 2.9 à 2.11.

2.2 Conceitos Iniciais para a Distribuição deSobrevivência

Chamaremos de T a variável aleatória que representaráo tempo de falha dentro do nosso estudo. Aqui a noção detempo é usada de maneira genérica. Ele pode serrealmentetempoou qualqueroutravariávelnãonegativa,realmentetempoou qualqueroutravariávelnãonegativa,desde que haja um número qualquer de falhas ou quebrasassociado a variável.Denotamos :sendo a distribuição de T e denotamos :

sendo a Função de Sobrevivência de T. Note que algunsautores definem F(t) e S(t) por Pr(T<=t) e Pr(T>t)respectivamente. Na prática isto não faz diferença para osresultados que se seguem quando T é uma variávelcontinua, este caso será considerado a partir de agora. Nósiremos assumir que T tem a função de densidade :

)Pr()( tTtF <=

)(1)Pr()( tFtTtS −=≥=


tal que a probabilidade da unidade falhar num curtoespaço de tempo [t , t+δt) é :

Considere a probalidade condicional do item falharnaquele instante [t , t+δt) e não ter falhado até o tempo t :

Podemos pensar como a probabilidade do itemiminentemente falhar em t. A função h(t) é dada por :

dt

tdS

dt

tdFtf

)()()( −==

ttftTt δδ )()Pr( ≅≤≤

)()(

)|Pr(tS

ttfttTt

δδ ≅≤≤

esta é a função de risco, de taxa de falha, ou hazard e éum indicador natural da propensão a falha após umaunidade de tempo ter transcorrido. A função de taxa defalha acumulada é dada por :

e concluímos que :

)()(

)(tS

tfth =

∫=t

duuhtH0

)()(

))(exp()( tHtS −= 2.1


Note que f, F, S, h e H são funções tais que oconhecimento de uma delas nos permite o cálculo detodas as outras.

Alguns casos típicos são discutidos aqui :• Se h(t)=λ é constante então H(t)=λt e S(t)=exp(-λt ) é adistribuição de sobrevivência exponencial com parâmetroλ. A densidade correspondente é f(t)=λ exp(-λt).• Se h(t) é uma função crescente de t , então T é dito teruma taxa de falha crescente (IFR). Isto á apropriadoquando a unidade medida tem relação com fadiga oudanos cumulativos.• Se h(t) é uma função decrescente de t, então T é dito terum taxa de falha decrescente (DFR). Isto pode ocorrer,por exemplo,quandoo processodiminui a quantidadepor exemplo,quandoo processodiminui a quantidadeproduzida ao longo do tempo diminuindo o risco de falha.Isto é comum em alguns ambientes de produção decomponentes eletrônicos.• Outro caso comum mencionado é o “bat-tub harzard”onde a função de taxa de falha é decrescente inicialmentee depois torna-se crescente. Isto costuma acontecer emlinha de produção onde os componentes iniciais tem umaqualidade melhor que os finais provocando este tipo deoscilação na taxa de falha.


2.3 A Distribuição Exponencial

Como mencionado na sessão 1.6, a distribuiçãoexponencial é o ponto natural de início para umadistribuição de sobrevivência. Relembrando temos que aDistribuição de Sobrevivência, hazard e função dedensidade tem a seguinte forma:

ondeλ é um parâmetro positivo, frequentemente chamado

),exp()(

)(

)exp()(

ttf

th

ttS

λλλ

λ

−==

−=2.2

de taxa, e onde t>0. Note também que a distribuiçãoexponencial tem média 1/λ e variância 1/λ2. A forma dadensidade é a mesma para todos osλ, e 1/λ age como umparâmetro de escala. Então, por exemplo, se o tempo desobrevivência, T, de um certo tipo de componente émedido em minutos e ele é distribuído exponencialmentecom taxa igualλ, então T*=T/60 medido em horas édistribuído exponencialmente com taxa 60λ. Uma outraformulação comum é termos a parametrizaçãoα= 1/λ nolugar de λ. A Figura 2.1 mostra duas densidades dadistribuição exponencial com diferentes taxas.


As funções hazard correspondentes são apresentadasna Figura 2.2. Nós iremos mostrar que a DistribuiçãoExponencial é um caso especial das Famílias deDistribuições Weibull e Gamma.


2.4 Distribuições Weibull e Gumbel

Uma variável aleatória Weibull (W. Weibull(1939,1951)) possui a seguinte função de sobrevivência:

para t>0 e ondeα e η são parâmetros positivos, sendoαum parâmetro de escala eη um parâmetro de forma. Noteque quandoη=1, obtemos uma Distribuição Exponencialcom parâmetroλ=1/ α.

)exp()(η

α

−= ttS 2.3


A função de falha (hazard) da Weibull é :

Então temos DFR paraη<1, constante paraη=1 e IFRparaη>1. Em particular, para 1<η<2 a função de falha seaproxima de uma função linear e paraη=2 a função élinear; paraη>2 a função cresce rapidamente acima deuma função linear. A função de taxa de falha (hazard) daWeibull para diferentes valores dos parâmetros é mostradana Figura 2.3.

A função de densidade da Weibull é

1)( −−= ηηηα tth

)exp()( 1η

ηη

αηα

−= −− tttf 2.4

para t > 0.

)exp()(α

ηα

−= ttf 2.4


A média e a variância são dadas por :

∫∞

−

−

−

−=Γ

+Γ=

+Γ=

0

1

1

1

)exp()(

,

)12(var

)1(

duuux

onde

iancia

media

x

ηαηα

2.5


veja , por exemplo, Abramowitz and Stegun (1972,CAPITULO 6). Um programa em fortran para calcular aequação 2.5 é dado em Griffiths and Hill (1985, pp. 243-6), que é baseado num programa anterior de Pike e Hill(1966). Quandoη é grande (>5), a média e a variância sãoaproximadamenteα e 1.64α2/ η respectivamente. A formada densidade depende deη. Na Figura 2.4 são mostradasalgumas funções de densidade da Weibull para diferentesvalores deη.


A Distribuição Weibull é provavelmente a mais utilizadadas distribuições em análise de sobrevivência. Umapossível explicação para isto se deve ao seucomportamento nos extremos da distribuição, àpossibilidade de variarmos o seu formato e em particulara possibilidade de utilizá-la como uma generalização daExponencial.

A Distribuição de Gumbel tem a seguinte função desobrevivência :

para , ondeµ é o parâmetro de locação eσ é oparâmetro de escala. Esta distribuição também começacom limite de distribuição mínimo, veja Galambos

])/)exp[(exp()( σµ−−= xxS 2.6

∞<<∞− x

com limite de distribuição mínimo, veja Galambos(1978), e tem uma taxa de falha exponencial crescente.Em alguns casos permite valores negativos comprobabilidades positivas. Mais comumente a distribuiçãode Gumbel é gerada através de Log(t) quando T tem umadistribuição Weibull. A relação entre os parâmetros daGumbel e da Weibull é a seguinte :

A função de densidade da Gumbel é a seguinte :

ησαµ

/1

)log(

==

)()/)exp(()( 1 xSxxf σµσ −= −


para , e tem a mesma forma para todos osparâmetros.

Note que a média e a variância da Gumbel sãoµ -γσ e(π2/6)σ2, respectivamente, ondeγ=0.5772 é a constante deEuler, e a distribuição é negativamente inclinada. Adensidade e a taxa de falha (harzard) para a distribuiçãode Gumbel comµ=0 e σ=1 é mostrada nas Figuras 2.5 e2.6 respectivamente.

∞<<∞− x


2.5 Distribuições Normal e Lognormal

A distribuição Normal é a distribuição mais comumenteutilizada em Estatística.

Em confiabilidade é geralmente usada como ummodelo para log T. A função de densidade dadistribuição lognormal é descrita pela equação abaixo:

σ)µ−−

πσ=

²

²(logexp

t²)t(f

22

1

{ }2/²exp)( σµ + =TE { }{ }12 −σσ + µ= ²)exp(²exp)T(Var

As funções de Sobrevivência e de Risco podem serescritas somente em termos de integrais.

Algumas densidades e funções de risco são plotadas nafigura 2.9 e 2.10. A função de Risco é crescente paravalores de t próximos de zero e eventualmentedecrescente quando .∞→t

)(a )(b


25.1)(1)(0)(2)(

1log49.2

= =.5=50.= ρρρρ dcba

emédiacomnormalõesdistribuiçparadensidadedeFunçõesFigura

)(a )(b

)(c )(d

25.1=ρ)d( 1=ρ)c(.50=ρ)b(520.= ρ)a(

e1médiacomnormallogõesdistribuiç4parariscodeFunções10.2figura

)(c )(d


2.6 Distribuições Gama e GamaGeneralizada

A distribuição Gama é descrita pela equaçãoabaixo:

Densidades são positivas (ver figura 2.11) mastendem para normal quando é grande. As funções deSobrevivênciae risco podem ser escritas somente em

2

1ρρ

σ

ρ=)T(Vare

λ

ρ=)T(E

)tλexp(t)p(Γ

λ=)t(f

Sobrevivênciae risco podem ser escritas somente emtermos de integrais.

A função de risco é decrescente para ,constante para (exponencial) e crescente para(ver figura 2.12).

A Gama é obtida como a distribuição doρ-ésimotempo de falha num processo de Poisson .

1<ρ

1=ρ 1>ρ


)(a )(b

)(b)(c )(d

5.0=ρ)d( 2.5=ρ)c(1.5=ρ)b(0.5= ρ)a(

e1médiacomgamaõesdistribuiç4paradensidadedeFunções11.2figura

)(a )(b


A distribuição gama generalizada é descrita pela equaçãoabaixo:

)(c )(d

5.0=ρ)d( 2.5=ρ)c(1.5=ρ)b(0.5= ρ)a(

e1médiacomgamaõesdistribuiç4parariscodeFunções12.2figura

abaixo:

A distribuição gama generalizada inclui os seguintescasos especiais:

{ }η1−ρη

)λ−)ρ(Γ)λ(ηλ= t(exp

t)t(f

1=ρ=η1=ρ

1 = η

:lExponencia)iii(

:Weibull)ii(

:Gama)i(


2.7 Distribuição Exponencial por Partes

Uma generalização da distribuição exponencial

Temos abaixo a função de Risco:

Vantagem: Pode-se aproximar qualquer função de

L,,),0( 21 IIemdepartiçãonaBaseado ∞

=∈

=∈=∈

=

1

],0[,

],0[,

],0[,

)( 222

11

mmm tIt

tIt

tIt

th

λ

λλ

M

Vantagem: Pode-se aproximar qualquer função deRisco desejada.

Desvantagem: Grande número de parâmetros (“não-paramétrica”)

{ }{ }

{ } iiii

m

ii

i

jjjji

iiii

iii

IttttTtT

kTE

kemittkonde

ttktf

miitttktS

,)(exp)/Pr(

)()(

1,...,2,)-(exp

)-(-exp)(

,...,1,,)-(-exp)(

11

11

11

2

11

1

1-

11-

1-

1-

λ

λλλ

λ

λλλ

=>>

+=

==

=

===

=

=∑


2.8 Censura

Observações incompletas frequentemente ocorrem nosestudos de sobrevivência e confiabilidade.

Nos testes de confiabilidade é comum aguardar atétodos os itens falharem.

Nos estudos de sobrevivência, pacientes abandonam otratamento ou continuam vivos depois do final dosestudos. Isso resulta em algumas observaçõesincompletas, ditas censuradas.

Tipos comuns de censura a direita:

Tipo I: Observações são acompanhadas até um tempo cfixado inicialmente.fixado inicialmente.Tipo II: Observações são acompanhadas até obter-se umnúmero pré-determinado de falhas.Tipo III: Aleatória à direita: Associado aos tempos defalha existem onde observa-se apenase

onde é o tempo de falha observado, eindependentes

iT s'Ci

>≤

=)(,0

)(,1

censuradoCT

censuradonãoCTX

ii

iii

iT s'Ci


2.9 Métodos dos Momentos para DadosSimples: SemCensura

Métodos informais (métodos formais serãoapresentados no próximo capítulo)

Baseado nos momentos e estimativas simplesSuponha que t1,...,tn sejam tempos de falha observados.

Exemplo 2.1: T – número de milhões de revoluções derolimã até a falha.

Dados: 17.88, 28.92, 33.00, 41.52, 42.12, 45.60, 48.40,51.84, 51.96, 54.12, 55.56, 67.80, 68.64, 68.64, 68.88,51.84, 51.96, 54.12, 55.56, 67.80, 68.64, 68.64, 68.88,84.12, 93.12, 98.64, 105.12, 105.84, 127.92, 128.04,173.40 (ordenados por conveniência)

Média amostral: = 72.22 e desvio padrão amostral (s.d.):st = 37.49Recai segundo §2.3 que a média e o s.d. coincidem nomodelo exponencial.Neste caso, /st se aproxima de 2, logo o modeloexponencial não é apropriado.Para ajustar Weibull e lognormal, é mais fácil trabalharcom xi = log ti.Novamente, média amostral: = 4.150 e s.d. amostral:sx = 0.534.

t

t

xRodrigues Fazenda 31

Estes cálculos valem para a médiaµ –γσ e s.d.πσ/√6da Gumbel, trazendo assim os momentos estimados =0.416 e = 4.390.

Em termos dos parâmetros da Weibull, temos: =exp( ) = 80.64

e = 1/ =2.40, diferente do 1.

De forma similar, os parâmetros estimados dalognormal são = 4.150 e = 0.534.

Outra aproximação é baseada na função desobrevivência empírica dada por

σ~

µ~

α~

µ~

η~ σ~

µ~ σ~

( )tS

( )n

YstdenúmerotS ti ≥= '__ˆ

é um estimador não paramétrico de S(t)

n possui distribuição binomial com média S(t) e ,segundo §1.15, um IC a 100%(1-α) para S(t) dado por

De forma similar, a função de taxa de falha acumuladaempírica é dada por

Onde o s.d. é dado por

( )tS( )tS

( ) ( ) ( )[ ] ( ) ( ) ( )[ ] 21

2

21

2

ˆ1ˆˆ,

ˆ1ˆˆ

−+

−−

n

tStSztS

n

tStSztS αα

( ) ( )tStH ˆlogˆ −=

( )( )

21

ˆ

ˆ1

−

tSn

tSRodrigues Fazenda 32

Os gráficos destas funções empíricas podem ser usadospara checar a adequação das hipóteses dos parâmetros.

Assuma os tempos ordenados t(1) < ... < t(n) .

• Salto de 1 / n tempo t(i).

• Realocado por 1- (i – 0.5) / n . (Outra forma possível)Modelo Weibull: S(t) = exp {- (t /α)η}

Log S(t) = - (t /α)η � log{-log S(t)}= η log t – η log α

Se o modelo Weibull é apropriado, então o gráfico de

é aproximadamente uma linha reta.

( )( )n

itS i

1--1ˆ =

( ) ( )( ){ }( )ii tSt ˆloglog,log −

Inicialmente os parâmetros estimados serão obtidos a partir de :

-η log α = intercepto

η = coeficiente angular

Modelo lognormal: ,Φ – funçãode distribuição de N(0,1).

Se o modelo lognormal for apropriado, então o gráficode (log t(i),Φ-1{(i-0.5)/n}) será uma linha reta.

( )

−Φ−=σ

µttS

log1

( ){ }σµ

σ−=−Φ − t

tSlog

11


Inicialmente os parâmetros estimados serão dados por :

-µ / σ= intercepto

1 / σ = coeficiente angular

Estes gráficos são dados nas figuras 2.13 e 2.14 doexemplo 2.1.Inicialmente os parâmetros estimados são (entreparênteses por momentos):

Weibull: η =2.3(2.4) eα = 77.3(80.6)

Lognormal:µ = 4.2(4.15) eσ = 0.56(0.53)

Diferentemente do baseado pelo método dosmomentos, o gráfico das probabilidades pode ser usadocom censura.

Elessãodefinidospor t < t , parar temposde falhaElessãodefinidospor t < t (r ), parar temposde falha(não censurados).

2.10 Estimador do Produto-Limite

O estimador do produto-limite (PL) ou de Kaplan-Meyer é um estimador não paramétrico da função desobrevivência.

Ele coincide com a função empírica de sobrevivênciaquando não há censura.


a1 < ... < ak – k tempos de falha distintos (a0 = 0)

d1, ...dk – número de falhas em cada tempo de falha (d0=0)

n1 < ... < nk – número de itens em risco em cada tempo defalha (nk = 0)

O estimador do PL é:

Esta é uma função escada começando do 1 para t = 0 ealterando-se a cada ak.

É como se a distribuição de falhas se concentrasse nospontos a1, ... , ak.

De acordo com a teoria assintótica , média e variânciade sãodadospor

( ){ }∏

<

−=

taj j

j

jn

dtS

:

1ˆ

( )tSde sãodadospor

S(t) e

H(t) pode ser estimado de forma similar por

De forma mais simples e intuitiva, podemos estimarH(t) usando

que é relacionado ao estimador .

Pode-se utilizar análise gráfica do estimador do PLpara avaliação da adequação de modelos Weibul e log-normal.

( )tS

( )tS ( ){ }∑

< −taj jjj

j

jdnn

d

:

( ) ( )tStH ˆlogˆ −=

( ) { }∑ <=

taj n

d

j j

jtH:

~

j

jn

djah =)(

~


Exemplo 2.3: Resistência de corda a uma certa tensão(em unidades codificadas).

Principais interesses:

• Qual a confiabilidade de uma corda após 53 unidadesde tensão ?

• O modelo de distribuição Weibull é apropriado ?

Da tabela 2.2 ,

e

Um IC de 95% para S(53) é dado por

(0.6849-1.69x0.0725, 0.6849+1.69x0.0725)=(0.54,0.83)

Fora 3 pontos isolados a figura 2.17 parece com umalinhareta.

6849.0)53(ˆ =S( ) ( )22 0725.00112.06849.0))53(ˆ( ==SVar

linhareta.

Investigação similar com modelo lognormal apresentaos mesmos resultados.

Tabela 2.1 Resistência de 48 cordas

Observações não censuradas36,3 41,7 43,9 49,9 50,1 50,8 51,9 52,1 52,3 52,352,4 52,6 52,7 53,1 53,6 53,6 53,9 53,9 54,1 54,654,8 54,8 55,1 55,4 55,9 56,0 56,1 56,5 56,9 57,157,1 57,3 57,7 57,8 58,1 58,9 59,0 59,1 59,6 60,460,7

Observações censuradas pela direita29,6 33,4 35,0 40,0 41,9 42,5


Tabela 2.2 Cálculo da amostral daExemplo 2.3

j aj nj dj (nj-dj)/nj S(aj+0) dj/(nj(nj-dj))0 48 0 1,0000 1,0000 0,00001 36,3 44 1 0,9773 0,9773 0,0005 2 41,7 42 1 0,9762 0,9540 0,0006 3 43,9 39 1 0,9744 0,9295 0,0007 4 49,9 38 1 0,9737 0,9051 0,0007 5 50,1 37 1 0,9730 0,8806 0,0008 6 50,8 36 1 0,9722 0,8562 0,0008 7 51,9 35 1 0,9714 0,8317 0,0008 8 52,1 34 1 0,9706 0,8072 0,0009 9 52,3 33 2 0,9394 0,7583 0,0020

10 52,4 31 1 0,9677 0,7338 0,0011 11 52,6 30 1 0,9667 0,7094 0,0011 12 52,7 29 1 0,9655 0,6849 0,0012 13 53,1 28 1 0,9643 0,6605 0,0013

A aproximação do erro padrão de S(53) , vem da equação (2.21) e tabela 2.2

W e ib u l l

3 .0 3 .5 4 .0 4 .5 5 .0

-3-1

1

W e ib u l l

lo g ( t )

log

(-lo

g(l1

))

3 .0 3 .5 4 .0 4 .5 5 .0

-20

12

L o g n o r m a l

lo g ( t )

l2


3. MÉTODOS ESTATÍSTICOS PARAAMOSTRAS SIMPLES

3.1 Introdução

Final do último capítulo: métodos estatísticos simples.Este capítulo: métodos mais formais, máxima

verossimilhança, inferência bayesiana dinâmica.

3.2 Estimação por Máxima Verossimilhança: Generalidades

Suponhaumaamostradetemposdevidas ntt ,...,1Suponhaumaamostradetemposdevidasde uma certa população.Todos os possuem densidade

, ondeCaso as observações não sejam censuradas então, a função de

verossimilhança é

para observações censuradas (a direita):A contribuição para a verossimilhança é a

probabilidade de sobrevivência após o tempo de censura.Separando os dados em conjuntos disjuntos: C - itens

censurados e U - para itens não censurados.

ntt ,...,1sti `

( )θ\tf ( )mθθθ ,...,1

=

( ) ( )∏=

=n

iitfL

1

\θθ

( ) ( ) ( )

= ∏∏

∈∈ Cii

Uii tStfL θθθ \\Rodrigues Fazenda 38

Para outras formas de censura existem outras expressões.

É mais conveniente trabalhar com Estimativa da máxima verossimilhança (EMV) de

Normalmente são obtidos resolvendo

Assumindo que q(p;θ) como o quantil 1-p de T, ouseja, Pr(T≥ q(p;θ) = S(q(p;θ)) = p

( ) ( )( )θθ Ll log=

( ) ( ).θlou ,θL omaximizand θ é θ

( )mj

l

j

,...,1 , 0 ==∂∂

θθ

( ) ( ) .φ de EMV o é θg=φ uma, a uma ação transformumafor θg=φ Se

( ) ( ). ˆp;q é p;q de EMV O θθ( )

( )m<kj,<1 ,

θθ∂

θl∂=J com J =V

∞→n quando ,V,θN~θ :aassintótic teoriaPela

jkkj

21-

. vcomo denotado e V matriz da k)(j, elemento o :Notação jk

( )

( )mkj,1 , J com J V

n quando ,V,N~ˆ :aassintótic teoriaPela

21- <<

∂∂−==

∞→

jkkj

l

θθθ

θθ


J é a informação observada da matrizEm particular,

O EMV possui muitas vantagens sobre todos os outrosmétodos clássicos de estimação:• Eleéuniversal;

. vcomo denotado e V matriz da k)(j, elemento o :Notação jk

( )

( ) ( )

( )11

2

m

1=j

m

1=kjk

kj

v θ∂

θg∂

à se-reduz φ de aassintótic a variância 1,=m um para ,particular Em

vθ∂

θg∂

∂

θg∂,φN~φ

:é oassintótic resultado o , θg = φ mação transforuma Para

∑∑

• Eleéuniversal;• Ele é invariante;• Ele possui boas propriedades assintóticas:

Consistência, normalidade assintótica e eficiência;• Distribuição assintótica é facilmente encontrada.

3.3 Máxima Verossimilhança (MV) estimação: ilustrações

( ) , tlog xsendo , observados são t,..., t vidade temposque Suponha iin1 =


( ) ( )

( )

( )0<

λ

r=

λ∂

λl∂

tλ

r=

λ∂

λl∂

tλλrlog = λl

:lExponencia ãoDistribuiç

.censurados sãor -n e censurados são não itensr Onde

22

2

n

1=ii

n

1=ii

∑

∑

λ aassintótic a variâncicom

t

r=λpor dado máximo de ponto um é 0=l′ para soluçãoA

2

n

1=ii∑

r aassintótic a variâncicom

( ) ( )

( )

( )σ

µxexpµx+µr+xσr=

σ∂

l∂

σ

µxexp+r=

µ∂

l∂

σ

µxexp

σµr

σ

x+σlogr=σ,µl

:Gumbel da formato noou

tα

1tlog1η+αlogηrηlogr=α,ηl

WeibullãoDistribuiç

in

1=ii

|U∈ii

n

1=i

i

n

1=i

i

U∈i

i

n

1=i

ηiη

U∈ii

∑∑

∑

∑∑

∑∑


µ depois e σ determinar para usados são numéricos Métodos

) µ envolve não ( 0=

σ

xexp

σ

xexpxi

σ+xr

1

σ

xexp

r

1logσ=µ

satisfaz EMV o Portanto,

∑

∑∑

∑

n

1=i

i

in

1=in

1=ii

n

1=i

i

µxµxl∂

σ

r=

µ∂

l∂

: são ) aassintótic a variância para usado ( EMV do derivada segundaA

n2

22

2

∑

σ

µxexp

σ

µxr=

σ∂

l∂

σ

µxexp

σ

µx=

σ∂µ∂

l∂

i

2n

1=i

i2

2

i

1=i

i

∑

∑

( )( )

( ) . σ,µN log que visto

ajustada, ãodistribuiç uma de através estimadosser podem quantis Os

0.56) & (0.534 0.522=σ e ) 4.2 & (4.150 4.150=µ :Normal-Log Modelo

0.435) & (0.416 0.476=σ e ) 4.35 & (4.390 4.405=µ : Gumbel Modelo

gráficas) sestimativa & (momentos

ML de sestimativa das osacompanhad 2.1 exemplo no dados Os :3.2 Exemplo

. λl deplot um mostra 3.1 Figura

0.137.=λs.e e 0.434=λ : lexponencia Modelo

23.05=t e 10=r 13=n possuem 2.2 exemplo do dados Os : 3.1 Exemplo ∑ i


Figura 3.1 Log-Verossimilhança para o tempo de vida de

Lambda

Log-

Ver

ossi

milh

ança

0.2 0.4 0.6 0.8 1.0 1.2

-26

-24

-22

-20

Figura 3.1 Log-Verossimilhança para o tempo de vida de um componente de um avião com distribuição Exponencial. A reta no gráfico foi feita para mostrar o intervalo de confiança 95% para lambda, baseado em W.

Note a grande diferença na calda

( ) ( )

) parenteses entre

estão (s.e quantis dos sestimativa as reúne abaixo 3.2 A tabela

1.0=1.282-Φ com σ1.282-µexp

é Normal-log ãodistribuiç uma de 0.1 quantil do EMV O

Quantil Weibull Log-Normalmediana 68.7(8.0) 63.4(6.9)

10% inferior 28.1(6.3) 32.5(4.8)


3.4 Intervalos de Confianças e Testes

θ = (θ1,...,θn) está divida dentro (θ(Α), θ(Β)) dasdimensões ma e mb.

Nós interessa testar a Hipótese H:θ(Α)= θ0(Α)

( , ) é um EMV de (θ(Α),θ(Β))(Α0) é um EMV deθ(Β) de H.Estão disponíveis dois procedimentos:1) Onde temos para H , que:

W(θ0(Α)) = 2{ι ( , ) − ι[θ0

(Α), (Α0)]} ∼ χ2( ma),é aproximadamente.

Grandes valores de W⇒ grandes diferenças emcomum com log – máxima verossimilhança⇒ grande

^)B(θ

^)A(θ

^)B(θ

^)B(θ

^)A(θ

suporte contra H.O teste da relação da MV rejeita H se W(θ0(Α)) > χα

2((ma)

ondeχα2( (ma) é 1 -α quantil deχα

2(ma).As regiões de confiança para são dadas por:

{ :W( θ(Α),χα2( ma)}

2) Seja VA a variância assintótica de .Então:W*( ) = ( − θ0

(Α))ΤVA-1( − θ0

(Α)),distribuição aproximadaχ2( ma) de H.

As regiões de confianças e os testes são obtidos a partirdas informações acima com W* substituindo W.

^)A(θ

^)A(θ

^)A(θ^

)A(θ^

)A(θ

^)A(θ


Em particular, paraθ(A) escalar, a região de confiançatorna-se:[ - z α/2VA

½, + z α/2VA½] ( um intervalo)

onde é usado o fato deχ2( 1) = [N(0,1)].Embora W e W* sejam assintoticamente equivalentes e

geralmente similares, preferimos W pela re-parametrização acima e por ser invariante.

Exemplo 3.1 (cont.): Temos o modelo exponencial,ondeθ(A) = λ, ma = 1 é W*(λ ) = 2{- 18.35 - 10logλ +23.05λ}O intervalo de confiança (IC = 95%) baseado em W* édado por: {λ:W(λ)£ 3.84}=[.22,.76] e como 3.84 =χ0.5

2(1). (Figura 3.1)

^)( Aθ

^)( Aθ

0.5

Então o intervalo de confiança (IC = 95%) baseado emW*é:[.434 – 1.96 x .137, .434+1.96 x .137] =[0.17,0.70](simétrico)

O intervalo de confiança (IC=95%) é dado por:[0.21,0.74] baseado em 2rλ/ ∼ χ2(2r).

Exemplo 3.2 (cont.): Desejamos testar a hipótese H deexponencialidade.

Usando o modelo Gumbel temos:θ(Α) = σ, θ(B) = µ, ma= 1 e H éσ = 1

W*(1) = 49.56 >> 3.84⇒ H0 rejeitado.Logo, W(1) = 15.50 confirma a rejeição de H como

esperado.Rodrigues Fazenda 45

O Exemplo 3.2 considera a hipótese de recursividade:onde a Exponencial é um caso especial da Gumbel.

São mais difíceis de considerar a hipótese de não –recursividade para o tratamento clássico estatístico.

3.5 Bondade do Ajuste

Enfoque formal: Encaixar o modelo dentro de umaclasse de modelos ( Exemplo 3.2) ou usar a formaexcelente.

Técnicas GráficasPlote o gráfico de QQ: Sejaµ e σ, parâmetros de

locaçãoe escala, (onde é estimadordo PL) e F é alocaçãoe escala, (onde é estimadordo PL) e F0 é afunção de distribuição paraµ = 0 e σ = 1.

O plote dos pontos [aj , F0−1(pj)] deveria ser linear.

Plote PP: junte os pontos (pj , F( aj , )). (esta linhadeveria ser y = x)

Pode ser usado fora do modelo de locação de escala.Plote SP: Para estabilizar a variabilidade de PP, plote a

transformação y = (2/π)sin-1x em ambos os eixos.As figuras 3.2 e 3.3 mostra os dados do exemplo 2.3

plotado em PP e SP.

^

θ


Goodness-of fit

Figura 3.2 Plote PP da Weibull para os dados da resistênciada madeira

0.00

0.20

0.40

0.60

0.80

1.00

0.00 0.20 0.40 0.60 0.80 1.00

Valores Observados

Val

ore

s E

sper

ado

s

da madeira

Figura 3.3 Plote SP da Weibull para os dados da resistênciada madeira

0.00

0.20

0.40

0.60

0.80

1.00

0.00 0.20 0.40 0.60 0.80 1.00

Valores Observados

Val

ore

s E

sper

ado

s


3.6 Elementos de Estatística Bayesiana

Incorpora informação subjetiva sobre o problema(experiência anterior).

É feita através da especificação de uma distribuição apriori P(θ).

Informação a priori vaga: a análise é guiada pelainformação dos dados.

Assuma, como antes, uma amostrat = ( t1, …, tn ) comdensidadef (t ;θ).

Isto é combinado com a priori e leva a

Fórmula);()(

)|( 1

tfptP

n

ii∏

==θθ

θde Bayes

Válido paraθ e t discreto e contínuo.P(θ | t) é a densidade aposteriori(dado os dadost).Como t é constante, A fórmula de Bayes pode ser

simplificada em

A constante removidaP(t) pode ser recuperada por

Estimativas a posteriori paraθ são obtidas através demedidas de locação deP(θ | t).

Exemplo: Considere os dados do exemplo 2.2 commodelo exponencial

)()|(

tptP =θ

)()()|( θθθ LptP ∝

θθθ dLptP ∫= )()()(


É conveniente atribuir a priori (distribuição gama).

Combinada com a verossimilhança de forma simples(priori conjugada).

Para especificar os valores de a e b assuma queacredita-se queλ está próximo de 0.5 e que é poucoprovável que ele seja menor que 0.2.

Então, tome a moda da priori igual a 0.5 e P(λ < 0.2)≈0.05⇒a = 3 e b = 4.

A posteriori é λ | t ~ Gamma(13, 27.05) com moda0.444 e média 0.481 (figura 6.1).

e)( -23.0510 λλλ =L

e)( b-1 λλλ −∝ ap

O desvio padrão a posteriori (priori) é 0.133 (0.175).

λ05.2712λ-23.0510λ-41-3 eλ∝ ) eλ)( eλ ( ∝ t)| λp(


Regiões de confiança são facilmente obtidas daposteriori.

Particularmente úteis são as regiões de maiordensidade a posteriori (HPD).

Por exemplo, o intervalo HPD de 95% para todoλ é[0.231, 0.758].

Interpretação da região HPD é simples diferentementedas regiões de confiança clássicas.

Inferência sobre funções paramétricas são obtidas demaneira similar.

Assuma interesse na confiabilidade num certo tempo t0

.Para o modelo exponencial isto éS ( t0 ; λ) = e-lt0, uma

funçãodeλ.0

funçãodeλ.A posteriori completa deS( t0 ; λ) pode ser obtida.Como exemplo, 0.01 =Pr (λ < 0.225 |t) = Pr (S ( t0 ;

λ) > e-0.225t0 | t) ⇒ a probabilidade a posteriori de queS(t0; λ) excedae-0.225t0 é 0.01.

Predição:assuma que se está interessado no tempo devida S de um novo item.

Inferência deve ser baseada na distribuição deS | ( t1,... ,tn ).

(S independente det dadoλ)Por exemplo, a densidade do tempo de vida de um

novo item é

d t)| p( )| p(s d t)| p( t), | p(s t)| p(s λλλλλλ ∫∫ ==


(a constante de proporcionalidade é 0.481).Informação a priori vaga é usualmente representada

por adequados valores pequenos dos parâmetros da prioriconjugada.

No exemplo, pequenos valores de a e b, como 0.5.O intervalo de 95% de confiança a priori é [0.001,

5.024]. (muito grande)A posteriori éGama(10.5, 23.55) com média 0.446 e

desvio padrão 0.138.O limite de uma distribuição a priori vaga é uma priori

não informativa.No exemplo,issoéobtidofazendoa,b → 0 => p(l) ∝l-

14

0

)05.27(1305.2712

0

)037.01()()( ) |p(s −∞

+−−∞

− +∝=∝ ∫∫ sdedee ss λλλλλλ λλλ

No exemplo,issoéobtidofazendoa,b → 0 => p(l) ∝l-1 (priori imprópria).

A posteriori éGama (10, 23.05) com média 0.434 edesvio padrão 0.137. (similar aos resultados da inferênciapor máxima verossimilhança).

A priori não informativa é um meio para obtenção daposteriori na ausência de informação a priori.


3.7 Outros Tópicos emInferência Bayesiana(§6.4 do livro texto)

Especificação de priorisNão é necessário que a especificação seja muito

precisa.Próxima seção: análise Bayesiana com especificação a

priori parcial.Verificação de inconsistências: análise pré-posteriori.

Análise conjugada é conveniente mas nem sempreapropriada.

Prioris não informativas devem ser usadas comcuidado: podelevaraabsurdos.cuidado: podelevaraabsurdos.

Parâmetros de distúrbioSuponha que os parâmetros são divididos em (ψ, φ)

ondeψ é de interesseφ é somente necessário (distúrbio) e pode ser eliminadofacilmente via

Densidade marginalExemplo: tempos de falha Weibull com S(t) = e-λtγ

A verossimilhança écom

φφ dtptp ∫ Ψ=Ψ )|,()|(

∏=

−=r

i

si

r etL1

)()(),( γλγλγγλ ∑=

=n

iits

1

)( γγ


Assuma que a priori é

Integrando com relação aλ temos

Entre chaves:p(t | γ) - verossimilhança marginal (ouintegrada) deγ.

Fator de BayesConsidere um teste Bayesiano de uma dada hipótese

H.Se H é uma região entãoPr(H | t) reflete a crença em

H a posteriori.Parao exemploem§3.6 e H : λ < 0.25, temosPr(H | t)

)]([

1

11 )()|,()(),( γλγλ λγγγλλγγλ sbr

ii

rarba etptpepp +−

=

−+−− ∏∝⇒=

∝ +=

+

∏

ra

r

ii

r

sb

t

ptp)]([

1)()|(γ

γ γ

γγ

Parao exemploem§3.6 e H : λ < 0.25, temosPr(H | t)= Pr( λ < 0.25 |t) = 0.2 => H é rejeitada

Similarmente, pode-se usar a razão de chances aposteriori dada por

A primeira razão do lado direito é a razão de chances apriori e o segunda é o fator de Bayes.

Isso representa a razão relativa de verossimilhançasentre as duas hipótesesH e .

Maiores valores do fator de Bayes⇒ maior apoio dosdados emH.

O fator de Bayes é útil quando deseja-se testarhipóteses nulas bilateraisH : θ = θ0

)|()Pr(

)|()Pr(

)|Pr(

)|Pr(

HtpH

HtpH

tH

tH =

H


Note quePr(H | t) será sempre 0.Há alguma controvérsia a respeito de quão adequado é

testar hipóteses nulas bilaterais.No exemplo, o fator de Bayes paraH : l = 0.25 é 0.81

⇒ a crença no valor 0.25 é reduzida pelos dados.

3.8 Modelos Bayesianos Dinâmicos

Baseado numa distribuição exponencial por partes paraos tempos de falência com suposição explícita de conexãoentre intervalos.

Distribuição E. P. (Exponencial por Partes): o risco éconstante nos intervalos.

É geralmente razoável assumir funções de riscoÉ geralmente razoável assumir funções de riscocontínuas⇒ algumas conexão entre valores de risco emintervalos sucessivos.

Forma matemática adotada para simplicidade: passeioaleatório na escala log

wi - termo de perturbação permitindo o aumento daincerteza como um movimento direto do intervalo.

Consequências deste modelo:(1) Preserva a posição:(2) Aumento da incerteza :

iiii1ii W=)w(Vare0=)w(Eondew+λlog=λlog


Um artifício útil para determinar valores para Wi´s:fatores de desconto controlam a quantidade de informação(medida de precisão) passando direto dos intervalos.

Fator de desconto é um número entre 0 e 1(geralmentefechado para 1).

1) se o desconto é fechado para 0 – nenhumainformação passa direto do intervalo. Estimação M.V.(máxima verossimilhança) com distribuição E.P.(exponencial por partes), estimador P.L. (produto-limiteou Kaplan Meyer).

2) se o desconto é 1 – toda informação passa diretodo intervalo

Parâmetros são os mesmos⇒ tempo de falênciaexponencial.exponencial.

Função de verossimilhança

De §2.7, a verossimilhança para um dado indivíduo é

onde Xi é o indicador da falência no intervalo Iibi é o tempo observado em Ii, para este indivíduo.

Para uma amostra de tamanhon, a verossimilhança é

}λbexp{λ iixii

∞

1=iΠ

i

∞

1=iL=L Π


onde

onde di é o número de indivíduos observados até a falha em Ii

ai é o tempo total observado em Ii para a amostraLi é a verossimilhança paraλi baseado nos eventos

observados em Ii dado Di-1, a informação do intervaloanterior (ver final de §2.7).De fato, o produto das verossimilhanças vai de i = 1 a i =N onde N é o indexador do último intervalo com tempode falência observado (censurado ou não censurado).

Análise Sequencial e Distribuições a Priori

}aλexp{λ=L iidii

i

Análise Sequencial e Distribuições a Priori

Assuma queA Priori assumida para

Consequentemente:

aumento da incerteza

Atualização da distribuição deλi feita direta da fórmulade Bayes.

)γ,α(G~Dλ 1i1i1ii

1≤c<0onde)γc,αc(GéDλ i1ii1ii1ii

)Dλ(E=)Dλ(E)1 1i1i1ii

)Dλ(Varc=)Dλ(Var)2 1i1i1

i1ii


e (análise conjugada)

As análises procedem do aumento de i para i+1 comoantes dito.

Inicia em i = 0 indo para i = N, o último intervalo cominformação de dados.{ci} controla a suavidade da função risco:ci→ 0: sem passagem de informaçãoci = 0: passagemtotaldeinformação

}λaexp{λ}λγcexp{λ∝ iidii1ii

1αci

i1ii

}λ)a+γc(exp{λ∝ ii1ii1d+αc

ii1ii

i1iiii L)Dλ(E∝)Dλ(p

)a+γc,d+αc(G~Dλ i1iii1iiii

ci = 0: passagemtotaldeinformação

Especificação dos ci´s

Dado queUsando o método Delta

Da relação entre sucessivosλ´s:

Usando novamente o método Delta dado

)γ,α(G~Dλ 1i1i1i1i

1i1i1i

1i

1i1i1i α

1=)Dλ(logVare)

γ

αlog(=)Dλ(logE

i1i

1i1i1i

1i1ii W+

α

1=)Dλ(logVare)

γ

αlog(=)Dλ(logE

)W+α1

()γ

α(=)Dλ(logVare)

γ

α(=)Dλ(E i

1i

2

1i

1i1ii

1i

1i1ii Rodrigues Fazenda 57

De onde um obtémGeralmente, Wi é selecionado proporcional ao

tamanho de Ii.Quanto maior o intervalo, mais informação é perdida.A proporcionalidade constante é W, a variância da

perturbação supera uma unidade de tempo.Fatores de desconto são associados diretamente com

W: tem que ser especificado só a primeira vez para ummodelo dado.

Inferência

Inferência é baseada na distribuição predita de um novo tempo de falência S baseado em DN, o (total) de

1i1ii )Wα+1(=c

novo tempo de falência S baseado em DN, o (total) de informação do dado.

Interesse particular: sobrevivência predita erisco predito

Estes são obtidos após a integração fora dosλ´s comrespeito a sua distribuição suavizada (ou filtrada)

Estas distribuições são obtidas via um algoritmorecursivo.

Modelo de Seleção

Um modelo é especificado pela escolha de :1) Priori para2) fator de desconto3) grade de intervalos

)Ds(S N

)Ds(h N

)Dλ(p Ni

01 Dλ


Modelos M1 e M2 podem ser comparados via seusfatores de Bayes

Cada verossimilhança marginal é obtida após aintegração fora dos parâmetros como segue

facilmente obtido.

)Mt(p

)Mt(p

2

1

∏N

1=i1ii )M,D}Iemobservadoseventos{(p=)Mt(p

i1ii

N

1=i1iii λd)M,Dλ(p)M,D,λ}Iemobservadoseventos{(p=∏∫

i1ii

N

1=ii λd)M,Dλ(pL=∏∫


Tempo de Falência do Sistema de TelecomunicaçãoPeríodo de Observação de 20 de Maio de 1985 a 31de Outubro de 1985

( Z = dias da instalação para a falência, cancelamento ou data de encerramento)c = censurado e u = não censurado

Z Z Z Z Z Z

164c 2u 45u 147c 139c 135c3u 155c 150c 101c 139c 135c164c 155c 150c 146c 139c 135c164c 155c 150c 1u 139c 1u163c 139u 149c 143c 138c 134c

163c 152c 149c 143c 40u 13u163c 152c 149c 143c 138c 134c163c 152c 149c 142c 138c 134c163c 152c 149c 10u 138c 134c163c 94u 149c 141c 138c 134c163c 94u 149c 141c 138c 134c

77u 151c 149c 141c 138c 133c162c 151c 149c 141c 138c 133c162c 151c 149c 34u 138c 133c73c 151c 115u 140c 138c 133c63u 151c 148c 140c 137c 133c

161c 151c 148c 140c 137c 64c160c 151c 147c 140c 137c 133c160c 151c 147c 140c 137c 133c67u 90c 147c 140c 137c 133c141c 151c 147c 140c 137c 133c

156c 151c 147c 54u 137c



3.9 O Estimador Atuarial

Utilizado em tabelas de mortalidade onde muitas vezesdados estão agrupados. Assume-se que a distribuição écontínua e divide-se o tempo em intervalos geralmenteiguais onde a taxa de falha é constante.

Ex: Em tabela de mortalidade, divide-se tempo (tempode vida de população) em intervalos 0 – 1, 1 – 5 (ou 0 –5), 5 – 10, 10 – 15, ... anos.

Raramente faz-se divisão ano a ano.

Suponha que a população temn indivíduos morrendonum ano com idadesy1, ..., yn. (Não há censura)

Se o indivíduo morre no intervaloi , sua contribuição àverossimilhança é:

A verossimilhança total é dada pelo produto dascontribuições individuais

( ) ( )11

1

1

−− −−−−−

=⋅Π iijjj tytt

i

ji ee λλλ


⇐ verossimilhança fatora em i

onde di = # de mortes no intervalo i

indivíduo k morre antes do intervalo i)

(indivíduo k morre no intervalo i)

(indivíduo k morre depois do intervalo i)

∑⋅ =

⋅−∞

=Π

n

kiki

i

i

xd

i

e 1

1

λλ

( ]

>−∈−

<=

−

−−

−

ikii

iikik

ik

ik

tytt

ttyty

ty

x

,

;,

,0

1

11

1

é o tempo total em risco no intervalo i

Por analogia, na estimação do modelo exponencial, o EMV de λi é:

Se todos os indivíduos morrem no final dos intervalos,

∑ ikx

∑=

∧

ik

i

x

diλ

( )∑ ⋅−= −k

iiiik rttx 1Rodrigues Fazenda 63

onder i = # de indivíduos observados no intervaloi.

Se também há censura, fórmulas não se alteram (apenasr i será diferente).

Vamos supor agora que a censura também está sujeitaa um mecanismo:

• Aleatório cuja taxa é constante ao longo dos mesmosintervalos

• Independente do mecanismo de falha (mortalidade)

Inferência acima não é alterada pela independência.Normalmente, em tabelas de mortalidade, dados são

fornecidosemformagrupada,istoé,sósãofornecidos:fornecidosemformagrupada,istoé,sósãofornecidos:

di = # de mortes no intervaloimi = # de censuras no intervaloi

Temos 3 grupos de indivíduos em cada intervalo:

i) r i – di – mi - sobrevivem ao intervaloiii) di - morrem no intervaloiiii) mi - são censurados no intervaloi

Vamos supor que taxa de censura éθi no intervaloi


Já vimos que verossimilhança fatora emverossimilhanças condicionais à história passada.

A contribuição dada à verossimilhança do intervaloide cada um dos 3 grupos acima dada sobrevivência até oinício do intervalo é dada por:

(i) Pr (Y > ti , C > ti | Y > ti-1 , C > ti-1)(ii) Pr (Y ≤ ti , C > Y | Y > ti-1 , C > ti-1)(iii) Pr (Y > C , C ≤ ti | Y > ti-1 , C > ti-1)

falha

(iii) (i)(iii) (i)ti

(ii)

ti-1 censurati

De fundo bi = t i – ti-1 temos:

(i) Pr(Y>ti | Y>ti-1).Pr(C>ti | C>ti-1) = exp{-biλi}exp{-biθi}

= exp{-bi(λi+θi)}Rodrigues Fazenda 65

(ii) ( ) ( )∫∫∞

−− >⋅>− y ic

t

t iY dcdytCcftYyfi

i11 ||

1

( ){ } ( ){ }∫∫∞

−− −−⋅−−=− y iii

t

t iii dcdytctyi

i11 expexp

1

θθλλ

( ){ } ( ){ }∫−

−− −−−−= i

i

t

t iiiii dytcty1

11 expexp θλλ

( )[ ]{ }iiiii

i b θλθλ

λ+−−

+= exp1

( )[ ]{ }iiiii

i b θλθλ

θ+−−

+exp1(iii) , por analogia a (ii)

Logo, a verossimilhança do intervaloi é dada por

( ) ( ) ( ) ( )[ ]iiiii

ii

ii

ii

iiiiiiii θ+λbexp1log{m+d+

θ+λθ

logm+θ+λ

λlogd+θ+λbmdr

Nenhum outro fator de verossimilhança depende deλi e θiEMV de λi e θi podem ser obtidos a partir da verossimilhançaacima dando

e( )

−−+

−=∧

i

iii

iii

ii r

mdr

mdb

dlogλ ( )

−−+

−=∧

i

iii

iii

ii r

mdr

mdb

mlogθ


Normalmente, é pequeno

Fazendo aproximação temos

Equivale a assumir que mortes e censuras se distribuemuniformemente nos intervalos.

A probabilidade de sobrevivência a um intervalo é

e pode ser estimada por

i

ii

r

md +

( )2

1log2x

xx +≅+

( )( ) |2 ii

i

iiii

ii rb

d

mdrb

d =+−

≅∧λ

( ) ii bii etYtS λ−− => 1|

|iiii rdb ee −− =

∧λ

No caso específico de uma tabela de mortalidade,n ébastante grande e a única informação éd1, d2, ..., ou seja,número de mortos em cada intervalo.

Podemos calcularr i pois , masxik nãosão fornecidos.

Os problemáticos são os indivíduos que morrem (ousão censurados) no intervalo.

Suposição: Dada a massa de indivíduos, é razoávelsupor que indivíduos morrem (ou são censurados)aleatoriamente no intervalo.

∑∞

=

=ij

ji dr


Logo, para esses indivíduos é tomada como

se não há censura

Como r i = r i+1 + di temos

ondeDi = { indivíduos que morrem no intervalo i}Si = {indivíduos que morrem após intervaloi}

22intervalo

no morrem queindivíduos de #

11

1 −−

− −⋅=

−+⋅

ii

iiii tt

dttt

( ) ( )1111

2 −+−− −⋅+

−−=+= ∑∑∑ iiiiii

iS

ikD

ikk

ik ttrttt

dxxxii

( ) ( )111

2 −+− −+−⋅= iii

iii ttr

ttd

A suposição acima é razoável:da é falha apenas nointervalo 0-1 ano onde a tendência à falha é nitidamentemaior perto de 0.

2

( )

+⋅−= +− 211i

iii

drtt

( )2

drtt= ii1ii

( ) |1 iii rtt ⋅−= −


Logo, λi é usualmente estimada por

A probabilidade de sobrevivência a um intervalo é

e pode ser estimada por

Muitas vezes é pequeno. Fazendo a expansão de

Taylor em torno de 0, obtemos

que é denotado por ou

( ) |1 iii

L

rtt

d

⋅− −

( ) ( )11| −−−

− => iii ttii etYtS λ

( ) |1 iiiii rdtt ee −−− =−

∧λ

|i

i

r

d

|1

|

i

ird

r

de ii −≈−

ip~ iq~1 −

e substituindo esses valores em

temos o estimador atuarial

Observe que o estimador atuarial difere do estimador PL pela troca de r i por quanto t = ti , i = 1, 2, ...

Nos outros pontos ele é contínuo e o PL é tipo escada1

0

( ) ( )∏=

−>=i

jjji tYtStS

11|

( ) ∏=

∧=

i

jji ptS

1

~


3.10 O Estimador Bayesiano

Quando não se assume nenhuma distribuição específicapara os tempos de falha a própria Fd F ou a f.s. S torna-seo parâmetro da distribuição.

Convenciona-se dizer que o problema é não paramétricopois a dimensão do parâmetro é infinita.Temos que construir priori sobre F(t) ( ou S(t) ), ∀ t

∈[0;∞)

Ferguson (1973): Sejaα uma medida finita R+, isto é,

α ( [a,b] ) = c > 0, α ( [a,b]∪A ) ≥ c e α ( [0,∞) ) < ∞

Ex: ∴ α ( [0,∞) ) = 1 < ∞

A distribuição P é um processo de Dirichlet se qualquerpartição B1, ..., Bk de R+, Pr(B1), ..., Pr(Bk) temdistribuição Dirichlet com parâmetro (α(B1), ...,α(Bk) )

( ) ∫−=

A

udueAα

( ) ( )kk D ααθθθ ,,~,, 11~

KK=Rodrigues Fazenda 70

onde α1 = α(A1), ..., αk = α(Ak)

e, portanto,

é uma amostra aleatória de um

processo de Dirichlet se:

Pr(X1∈A1, ..., Xn∈An | P(A1), ..., P(An) ) =

Pode se mostrar que

∏=

−∝k

ii

ip1

~

1)( αθθ

( )nXXX ,,1~

K=

( )∏=

n

iiAP

1

( ) ( )=∈Pr 1α AAXPode se mostrar que

O processo de Dirichlet funciona como a priori para a distribuição amostral. Assim, por exemplo,

Logo tem densidade

e

α corresponde ao valor esperado a priori para F.

( ) ( )[ )( )∞

=∈;0

Pr 111 α

α AAX

( ) [ ]( ) ( )21;~;0Pr ααDttF =

( ) [ ] [ ] 11 21 )(1)()( −− −∝ αα tFtFtFf

[ ] [ ]( )[ )( )∞

=+

=;0;0

)(21

1

αα

ααα t

tFE


Se não há censura na amostra, a distribuição éconjugada e a posteriori de P também é um processo deDirichlet com parâmetroβ onde:

, se α é contínuo, β não é mais.

Se há censura, a distribuição não é conjugada e a forma da posteriori complica. Pode-se obter a esperança a posteriori de F(t) ou S(t), ∀ t > 0

Usando essa esperança como estimador temos

( ) ( ) ( )∑=

=+=n

ii tyItt

1

αβ

)(tS∧

( )( ) ( ) ( )( ) ( )( )( ) ( )∏

∧ +∞⋅++∞=l

ii yrytrttS

;0;)(

αα

yl < t ≤ yl+1c = k, ..., mondey1, ..., yk são as observações não-censuradasyk+1, ..., ym são os diferentes valores das observações censuradasck+1, ..., cm são o número de observações censuradas = yk+1, ..., ymr(t) = # de observações à vista em t

Se não há censura,

yi-1 ≤ t < yi , i = 1, ..., n

( )( ) ( )[ )( )

( )( ) ( )( )( ) ( )∏

+= −+∞+∞⋅

+∞++∞=

kj jjj

ii

cyry

yry

n

trttS

1 ;

;

;0

0;)(

αα

αα

( )( )[ )( ) n

jnttS

+∞+−+∞=

∧

;01;

)(α

α


No caso geral, se α→0, → estimador PL

Quando se α > 0, comporta-se da seguinte forma

)(tS∧

)(tS∧

A(n)

Abordagem completamente não paramétrica sugerida por Hill (1968) baseada na hipótese A(n):

Sejam X1, ..., Xn+1 permutáveis com distribuição P.

Suponha que observa-se X1 = x1, ..., Xn = xn ondex1 < x2< ... < xn (possível pela permutabilidade)

Sejam I(0), ..., I(n) intervalos dados por I(i) = (xi, xi+1) onde x0 = 0 e xn+1 = ∞


Então a distribuição preditiva de Xn+1 dado X1 = x1, ...,

Xn = xn dá

Eventuais empates nos xi’s podem ser separados acrescentando ε pequeno a um deles

No caso de censuraa situaçãocomplica pois não

( )1

1)(Pr 1 +

=∈+ niIX n

No caso de censuraa situaçãocomplica pois nãoconhecemos todos osxi’s. Resolve-se também de umaforma não paramétrica.

Suponha

x4 ex5 são tempos de censura.

O objectivo é calcular as probabilidades preditivas paraI(0), I(1), I(2), I(3). Isso é feito considerando todos ospossíveis valores (não censurados) dex4 ex5.


Supondo que x4 ∈ I(1) e x5 ∈ I(2) temos sob a hipótese A(5) que a probabilidade

de I(0) é , de I(1) é , de I(2) é e de I(3) é

Se supomos, por exemplo,

61

31

61

61 =+ 3

16

16

1 =+

61


Implícito nos cálculos acima:

Falha dos censurados poderia ocorrer em qualquerponto no intervalo de censura, ou seja, o ponto de censuraé trazido de volta até o início do intervalo.

Ex.: x4 poderia ocorrer em qualquer ponto de I(1)O cálculo é portanto uma aproximação (que fornece cotasuperior).DefinaZ1, ..., ZN – tempos de falha a ser observados

Z – tempos de falha a ser previston < N – número de falhas observadas

– valoresobservadosdefalha( )xxx ,,K= – valoresobservadosdefalha

yn+1, ..., yN – valores observados de censura

IEC = {Zj > yj, j = n+1, ..., N} Info. Exata de CensuraIPC = { Zj > U j, j = n+1, ..., N} Info. Parcial de Censura

onde Ui é o maiorxi | xi < yi

Ex.: No exemplo, temos x1 < y4 < x2 < y5 < x3 ∴ U4 = x1e U5 = x2

Ao invés de calcularmos Qi = Pr( Z ∈ I(i) | IEC, ), i = 0, 1, ..., n calcularemos

( )nxxx ,,1~

K=


Pi = Pr( Z ∈ I(i) | IPC,

Como IEC ⊂ IPC, Qi = Pr( Z ∈ I(i) | IEC, IPC, ) =

Como, dado IPC, pouca informação adicional sobreIEC é fonecida por [Z∈ I(i)] principalmente em amostrasmoderadasougrandes,Q ≈ P.

~x

) x IPC, |Pr(IEC

) x IPC, |IEC I(i), ZPr(

~

~∈

=

iP⋅∈

=) x IPC, |Pr(IEC

) I(i), Z,x IPC, |IEC Pr(

~

~

moderadasougrandes,Qi ≈ Pi.

Para o cálculo dos Pi, define-se ci - # de observaçõescensuradas em I(i)

, i = 0, 1, ..., n

λi = 1 / (N - (i - 1) - Ci)

Então P0 = λ0 e , i = 0, 1, ..., n-1

∑=

=i

kki cC

0

( )∏=

++ −⋅=i

jjiiP

111 1 λλ


Prova: (i = 0): Se c0 = 0, C0 = 0 e A(N)⇒ P0 = 1 / (N+1)

Se c0 ≠ 0, sob IPC, essas c0 observações são colocadasno início de I(0) e portanto não trazem nenhuma info.

Podemos considerar apenas N-c0 = N-C0 observaçõesrestantes.

A(N-C0) ⇒ P0 = 1 / (N-C0+1)

Sob IPC, c1 observações são colocadas no início de I(1) e para o cálculo da parcela acima (condicional a Z > x1) não trazem nenhuma info. assim como x1.

Podemos então sob A(N – 1 – c0 – c1) = A(N – 1 – C1)

obter Pr( Z ∈ I(1) | Z > x1 , IPC, ) = λ1~x


Logo,

O mesmo raciocínio pode ser seguido e usando indução mostra-se o resultado.

A substituição de IEC por IPC faz com que

λi = Pr( Z ∈ I(i) | Z > xi , IPC, ) produzam cotassuperiores para Pr( Z∈ I(i) | Z > xi , IEC, )

Observe que

( )011 1 λλ −⋅=P

~x

~x

Cotas inferiores podem ser obtidas se censuras sãotrazidas para cima.

IPCI = {Z j > Uj , j = n+1 , ..., N} Info. Parcial deCensura Inferior onde Uj é o menorxi | xi > yj

Nesse caso Q0 ≥ Pr( Z ∈ I(0) | IPCI, ) ,por A(N) 1

1

+=

N~x


Assim como

Podemos definir Pr( Z ∈ I(i) | Z > xi , IPCI, )

Por analogia com temos que , i = 0, 1, ..., n

com C-1 = 0

Definindo então Pi = , = Pr( Z ∈ I(i) | IPCI, ) temos

( ) ( )i

1=j

sj

s1+i

i

1=jj1+i1+i λ1λ=λ1λ=P

=Iiλ ~

x

siλ

1)1(1

−−−−=

i

Ii CiN

λ

siP I

iP~x

Além disso,

= Pr( Z ∈ I(1) | Z > x1 , IPCI, ) ≤ Pr( Z ∈ I(1) | Z > x1 , IEC, ) ≤

Multiplicando as duas inequações tem-se

I1λ

~x

~x s

1λ


Em geral


4. MODELOS DE REGRESSÃO PARADADOS DE CONFIABILIDADE

4.1 Introdução

Até agora, itens pertencem a mesma população.Às vezes, outras variáveis afetam os tempos de falha.Exemplos: tensão, pressão, temperatura (confiabilida-

de), idade, tratamento, sexo (análise de sobrevivência).Estas variáveis são chamadas covariáveis, e devem ser

incorporadas ao modelo.Elas podem ser contínuas (tensão, pressão,

temperatura,idade)oudiscretas(tratamento,sexo).temperatura,idade)oudiscretas(tratamento,sexo).De agora em diante, analisaremos dados de falha com

covariáveis (regressão).Em confiabilidade, a maioria dos modelos são

baseados na distribuição Weibull.Outra opção é a lognormal: após a transformação log

nos tempos de falha, podemos fazer uso da teoria normale da regressão padrão.

Em geral, prefere-se a distribuição Weibull devido afacilidade do seu uso com dados censurados e devido aforma da função taxa de falha.

Propósito do estudo: determinar o quantoT é afetadoporx (covariáveis).

Exemplo:x – tensão;T pode decrescer comx.§ 4.2 a 4.6 descrevem diferentes modelos.


§ 4.7 a 4.9 lidamos com modelos baseados na distribuiçãoWeibull.

Capítulo 5 lida com modelos de taxa de falhaproporcional.

Capítulo 6 lida com modelos Bayesianos dinâmicos.

4.2 Modelos de Tempo de Vida Acelerado(ALM)

Suponhamos tempos de falha sujeitos a uma carga.ALM : tempo de falha é o produto de uma função da

carga e do tempo de falha padrão.Tempode falha padrão: tempode falha paraum nívelTempode falha padrão: tempode falha paraum nível

padrão de cargaPr (T ≥ t | x) = S (t; x) = S0 (t ψx )

aondeS0 é a sobrevivência básica (padrão) eψx é umafunção positiva dex.

Quandox está no nível básico:ψx = 1. Para outrosvalores dex, tempo de falha é acelerado (multiplicado)porψx.

Exemplo:S0 é exponencial unidade (S0 (t) = e − t ) eψx = x β

⇒ S(t; x) = S0 (t x β ) = exp{ −t x β }

No ALM: Tyx tem distribuição básica oulog T = − log

yx + log W, aondeW ~ P(não depende dex).Rodrigues Fazenda 83

Uma especificação comum élog yx = xTb (similar aosmodelos lineares normais).

A função taxa de falhah (t; x) é

A adequação do ALMpode ser preliminarmenteavaliada pelos seguintes gráficos:

1) Sendolog T = − log yx + “termo de erro”, o gráficode log T versusx deve fornecer uma indicação da formadey.

2) Se os dados podem ser agrupados emk gruposhomogêneos, os gráficos das estimativas das funçõesde sobrevivência versus log t devem ser cópias

( ) ( ) ( )=

∂∂

−=∂

∂−=t

tS

t

xtSxth xψ0log;log

; ψx h0 (t ψx )

de sobrevivência versus log t devem ser cópiashorizontais deslocadas deS0 quando

Sj (t) = S0 (t ψx) = S0 { exp(log t + log ψj )} = s0 (log t + log ψj )

3) Também no ALMcom k grupos, os quantis são

proporcionais. Isto é, seqj (p) é o p-ésimo quantil do

grupo j entãoqj (p) ψj = qk (p) ψk porqueS0 (qj (p) ψj) =

Sj (qj (p)) = p, ∀j.

Os quantis 0.1, 0.2, ..., 0.9 para todos os grupos podem

ser estimados de .Rodrigues Fazenda 84

Podemos fazer o gráfico dos quantis estimados para o

grupo j, j = 2, ..., k versus os do grupo 1. Os gráficos

devem ser aproximadamente lineares com inclinaçãoψ1

/ψj .Equivalentemente, os gráficos dos log-quantis devem

ser linhas paralelas.Os gráficos também poderão ser feitos versus a média

dos quantis (ao invés de versus os quantis do grupo 1).

4.3 Modelos de Taxas de Falha Proporcionais(PHM)

PHM : a taxade falha é acelerada,isto é, h (t; x) = hPHM : a taxade falha é acelerada,isto é, h (t; x) = h0

(t) ψx aondeh0 é a função taxa de falha básica (padrão) eψx é uma função positiva dex.

Quandox está no nível básico:ψx = 1.

Este nome vem do fato (não depende det).

Uma especificação comum élog yx = xTb. Assim,h (t;x) = h0 (t) exp{ xTb}.

A função de sobrevivência é

( )( )

2

1

2

1

;

;

x

x

xth

xth

ψψ

=

S (t; x) = ( ){ }∫−t

dtxth0

;exp = ( ){ }∫−t

x dtth0 0exp ψ = ( )[ ] xtS ψ

0


SeS(t; x) é um ALM e um PHMentão

S(t; x) = Sa 0 (t y ax) =

e a única solução possível é a sobrevivência da Weibull

Sa 0 (t) = Sp 0 (t) = exp {−t h } com .

Gráficos preliminares podem novamente serem feitos

depois da separação dos dados emk grupos.

Assim, temos que

−log Sj (t) = −y j log S0 (t) ou log { −log Sj (t)} = log y j +

log { −log S0 (t)}

e os gráficos de devem ser múltiplos mutuamente na

escalalog ou devem ser cópias verticais deslocadas de

( ){ } pxtSpψ

0

ηψψ axpx =

jS

cada um na escalalog-log.

4.4 Modelos de Razão de ChancesProporcionais (POM)

POM: a função taxa de falha satisfaz

Após diferenciar em relação a t,

A razão taxa de falha éψx quandot = 0 e tende a 1quandot → ∞.

( )( )

( )( )tS

tS

xtS

xtSx

0

01

;

;1 −=− ψ

( )( )

( )( )

( )( )tS

xtS

tS

xtS

th

xthx

000 1

;1;;

−−==ψ


⇒ ocorre diminuição do efeito dex na taxa de falha a

medida que o tempo cresce.

Avaliações preliminares podem ser feitas baseadas nos

gráficos .

4.5 Generalizações

A relação PHMlog { −log Sj (t)} = log ψ j + log { −log S0 (t)}

pode ser generalizada comoφ 1 { −log Sj (t)} = g 1 (ψ j ) + φ 1 { −log S0 (t)}.

Isto inclui o POMcomo o caso especial:φ 1 (u) = log{( 1 − u) / u} .

{ } jj SS ˆˆ1−

{( 1 − u) / u} .A relação ALM log q j = log q0 − log (ψ j / ψ 0) pode

ser generalizada comoφ 2 (q j ) = φ 2 (q0) −g2 (ψ j / ψ 0).Outras generalizações são:1) ALM generalizado aondey x também depende det.2) PHM generalizado aonde aondes (l ) é a família

Box-Cox de transformações ems.

3) Modelo de deslocamento no tempo:h (t; x) = h0 (t +

yx) – taxa de falha de ação atrasada.

4) Modelo de taxa de falha polinomial:h (t; x) = y 0x +

y 1x t + ... +y qx t q.

5) Modelo de taxa de falha por partes:h (t; x) = h i (t;

x), parat ∈ I i, i = 1, ...,N.Rodrigues Fazenda 87

6) Covariáveis dependentes do tempo:x pode dependerdo tempo em alguns casos.

4.6 Modelos Bayesianos dinâmicos (DMB) (cf Gamerman, 1991)

DMB combina modelos de taxas de falhaproporcionais (PHM) com distribuição exponencial porpartes (similar ao modelo de taxa de falha por partes).

Do PHM, temos queh (t; x) = h0 (t) y x.Assuma queh0 (t) = exp{ b i, 0}, t ∈ I i.Também, generalizey x paray i, x, para que dependa de

t.t.Uma especificação comum éy i, x = exp{ xTb i}.O modelo pode ser escrito comoh (t; x) = exp{ xTb i},

aondeb i agora inclui uma primeira componenteb i, 0 e xT

agora inclui uma primeira componente 1.O modelo é completado com a relação entre os

parâmetros em intervalos sucessivos, como em § 3.8.b i = G i b i−1 + w i aonde E(w i) = 0 e Var (w i) = Wi.A matriz evoluçãoG i fixa a parte determinística da

evolução e o termo de errow i controla o aumento naincerteza a medida que o tempo passa.

Exemplos:1) G i = I, é a matriz identidade (passeio aleatório

simples); Rodrigues Fazenda 88

2) G i = diag (1, G i 1) aonde

(modelo de crescimento generalizado).A variância pode novamente ser especificada através

dos fatores de desconto.Pode ter uma para cada parâmetro mas o mais comum

é ter uma para o parâmetro de baseb i, 0 e uma para oscoeficientes de regressão.

Geralmente, esta última é próxima a 1 (se 1, o modelose torna PHM).

Se não há presença de covariáveis, o modelo se torna omesmo de § 3.8.

=

10

)(11

ii

ItamanhoG


4.7 Modelos Baseados na DistribuiçãoWeibull

O Modelo de Regressão Weibull pode ser escrito como:

( )

−=

η

α x

txtS exp; onde βα T

x x=log

}exp{)(0ηttS −=

⇓1

Pode, também, ser escrito como:

( ) )(; 0 axtSxtS ψ= ondex

ax αψ 1= ALM

[ ] pxtS ψ)(0= onde ηαψ

xpx

1= PHM

WT x ηα 1

loglog +=


Sabemos que W ~ Gumbel e ainda queη é somente um fatorescala, mas também pode ser dependente de x.

Estimação e Testes

Regressão Linear Simples:

Estimador de Mínimos Quadrados pode ser usado.Procedimento similar pode ser usado para regressões múltiplas.Mais formalmente, o Estimador de Máxima Verossimilhançapode ser calculado como no capítulo 3: parâmetros sãoestimados, variâncias assintóticas obtidas e testes de hipótesescalculadosvia TestedaRazãodeMáximaVerossimilhança.

onde: ( )w

w

xx

eeTwT

wW −=

≥=

≥

=≥ expPrlogPr)Pr( ηα

ηα

xx 10log ββα +=

calculadosvia TestedaRazãodeMáximaVerossimilhança.

Esses cálculos podem ser manuseados por GLIM.

Plot de ResíduosPPplot pode ser usado como em 3.5

Os dados transformados convergem para uma amostra de variáveis aleatórias U(0,1).Os ui’s são chamados resíduos generalizados.Também, zi = log(-log ui) tem distribuição Gumbel.Um plot de probabilidade pode ser construído como segue:

( )φ,; iii xtSu =


• Ordene os zi’s em zi:n, ..., zn:n

• Faça o gráfico zi:n contra log{- log(1-pi:n)} onde pi:n é:

Dados não censurados

• Checar se o gráfico é da forma y=x

Resíduos generalizados são também proveitosos quandoplotados contra covariáveis (ver 4.8)

n

iou

n

5.01 −

∏∈≤ +−

+−+−Jjij jn

jn

n

n

, 5.1

5.05.01 Dados censurados


4.8 Um Exemplo: Resistências de Fibras de Carbono e PacotesDados das tabelas 4.1 e 4.2 são analisados.

2,247 2,640 2,842 2,908 3,099 3,126 3,245 3,328 3,355 3,383 3,572 3,581 3,6813,726 3,727 3,728 3,783 3,785 3,786 3,896 3,912 3,964 4,050 4,063 4,082 4,3114,118 4,141 4,216 4,251 4,262 4,326 4,402 4,457 4,466 4,519 4,542 4,555 4,6844,632 4,634 4,636 4,678 4,698 4,738 4,832 4,924 5,043 5,099 5,134 5,359 5,4735,571 5,684 5,721 5,998 6,060

1,901 2,132 2,203 2,228 2,257 2,350 2,361 2,396 2,397 2,445 2,454 2,454 2,4842,518 2,522 2,525 2,532 2,575 2,614 2,616 2,618 2,624 2,659 2,675 2,738 2,740

Tabela 4.1(a) Length 1mm

(b) Length 10mm

2,856 2,917 2,928 2,937 2,937 2,977 2,996 3,030 3,125 3,139 3,145 3,220 3,2233,235 3,243 3,264 3,272 3,294 3,332 3,346 3,377 3,408 3,435 3,493 3,501 3,5373,554 3,562 3,628 3,852 3,871 3,886 3,971 4,024 4,027 4,225 4,395 5,020

1,312 1,314 1,479 1,552 1,700 1,803 1,861 1,865 1,944 1,958 1,966 1,997 2,0862,021 2,027 2,055 2,063 2,098 2,140 2,179 2,224 2,240 2,253 2,270 2,272 2,2742,301 2,301 2,339 2,359 2,382 2,382 2,426 2,434 2,435 2,478 2,490 2,511 2,5142,535 2,554 2,566 2,570 2,586 2,629 2,633 2,642 2,648 2,684 2,697 2,726 2,7702,773 2,800 2,809 2,818 2,821 2,848 2,880 2,954 3,012 3,067 3,084 3,090 3,0963,128 3,233 3,433 3,585 3,585

1,339 1,434 1,549 1,574 1,589 1,613 1,746 1,753 1,764 1,807 1,812 1,84 1,8521,852 1,862 1,864 1,931 1,952 1,974 2,019 2,051 2,055 2,058 2,088 2,125 2,1622,171 2,172 2,18 2,194 2,211 2,27 2,272 2,28 2,299 2,308 2,335 2,349 2,3562,386 2,39 2,41 2,43 2,431 2,458 2,471 2,497 2,514 2,558 2,577 2,593 2,6012,604 2,62 2,633 2,67 2,682 2,699 2,705 2,735 2,785 2,785 3,02 3,042 3,1163,174

(d) Length 50mm

(c) Length 20mm


2,526 2,546 2,628 2,628 2,669 2,669 2,710 2,731 2,731 2,731 2,752 2,752 2,7932,834 2,834 2,854 2,875 2,875 2,895 2,916 2,916 2,957 2,977 2,998 3,060 3,0603,060 3,080

2,485 2,526 2,546 2,546 2,567 2,628 2,649 2,669 2,710 2,731 2,752 2,772 2,7932,793 2,813 2,813 2,854 2,854 2,854 2,895 2,916 2,936 2,936 2,957 2,957 3,0183,039 3,039 3,039 3,080

2,110 2,260 2,340 2,440 2,510 2,510 2,570 2,570 2,610 2,610 2,610 2,650 2,670(c) Length 150mm

Tabela 4.2(a) Length 20mm

(b) Length 50mm

2,110 2,260 2,340 2,440 2,510 2,510 2,570 2,570 2,610 2,610 2,610 2,650 2,6702,710 2,710 2,710 2,750 2,750 2,750 2,750 2,770 2,770 2,790 2,830 2,830 2,8302,870 2,870 2,900 2,900 2,920 2,940

1,889 2,115 2,177 2,259 2,279 2,320 2,341 2,341 2,382 2,382 2,402 2,443 2,4642,485 2,505 2,505 2,526 2,587 2,608 2,649 2,669 2,690 2,690 2,710 2,751 2,7512,854 2,854 2,875

(d) Length 300mm


Sob ALM, a distribuição de logT:

1) Variam em locação de um pro outro.2) Tem a mesma variância.

Variância amostral dos dados da tabela 4.1são:0.042, 0.040, 0.045 e 0.037 => estaticamente equivalentes

figura 4.1 é o gráfico 2 da seção 4.2: variação horizontal.

Figuras 4.2e 4.3: linear (particularmente na figura 4.3).Gráficos de contralog t deve ser:

•Cópiasverticaisdamesmacurva,sobPHM;

))ˆlog(log( jS−

•Cópiasverticaisdamesmacurva,sobPHM;•Linhas paralelas, sob modelo de regressão Weibull (verseção 4.7)

figura 4.4 mostra gráfico de contralog tpara dados da tabela 4.1com linhasWeibull (ajuste por Máxima Verossimilhançaseparadamente em cada amostra).

Figura 4.6 mostra o mesmo gráfico para os dados databela 4.2

Sob POM, devem ser cópias verticais uma da outra (verfigura 4.5)

))ˆlog(log( jS−


Agora, nesta Seção, vamos nos concentrar em WeakestLink Hypotesis:

onde Sr é a função de sobrevivência para fibra decomprimento r.

WLH pode ser implementada nos seguintes modelosWeibull:

LL tStSWLH )]([)(: 1=

−=

−=−

ηη

LL t

t

t

tLtSM expexp)(

10 onde η

1

1

−= LttL

10, ≤≤=−

εηε

Ltt

M1 – como em M0 mas com

M2 – como em M1 mas com

(modelo log-linear no parâmetro escala da distribuiçãoGumbel);

M3 – distribuição Weibull separada para cadacomprimento.

10,1 ≤≤= εηLttL

Llogloglog 1 τηη +=


A função de verossimilhança pode ser obtida para cada modelo e na maioria dos casos, maximizá-la não é difícil.

Os valores maximizados da função log-verossimilhaça são:

l0 = -229.1 (2 parâmetros)l1 = -227.7 (3 parâmetros)l2 = -227.6 (4 parâmetros)l3 = -220.1 (8 parâmetros)

Teste de hipótese H1: ε=1 é calculado via estatística de razão de verossimilhança:

é aceita.

(O Estimador de Máxima Verossimilhança sob M1 é 0.90)

Dados da tabela 4.2: valores maximizados da log-verossimilhança são:

84.3)23(8.24.12)(2 205.0011 =−<==−= χxllW

1H⇒


l0 = 21.8 (2 parâmetros)l1 = 29.6 (3 parâmetros)l2 = 31.5 (4 parâmetros)l3 = 35.8 (8 parâmetros)

Teste de H1 agora tem W1 = 15.6 > 3.84 => M0 rejeitado

O Estimador de Máxima Verossimilhança de ε sob M1 é 0.58

WLH não necessariamente implica modelo Weibull

PHM também é possível.





4.9 Análise Bayesiana de Dados de ConfiançaComo mencionado anteriormente, as análises de grande

amostra de máximo verossimilhança oferecem resultadossimilares à análise bayesiana com vaga informação a priori.

Exemplo: A estimativa de máximo verossimilhança 0,58 deξpara os dados da tabela 4.2 é uma aproximação da média daposteriori deξ.

Os testes de níveis de significância não possuemcomplemento bayesiano.

O teste de hipótese bayesiano pode ser feito da seguintemaneira:

1) Constrói-se a região de máxima densidadea posteriori1) Constrói-se a região de máxima densidadea posteriori100(1-α)% do parâmetro de interesse (ξ, no exemploanteriormente citado);

2) Observa-se se esta região contém o valor do parâmetroespecificado porH (1, no exemplo anterior);

3) Aceita-seH se este for o caso ou rejeita-seH, caso contrário.

Também, para grandes amostras, a razão de máximaposteriori pode ser aproximada pela razão de máximaverossimilhança, se a priori for vaga.

Exemplo: Considere os dados da tabela 6.1 do modeloWeibull.

onde W~ GumbelWatemperaturt 2/11)2.273(log −− +++= φβαRodrigues Fazenda 102


5. TAXA DE FALHA PROPORCIONAL

5.1 Introdução

Baseada no documento de Cox (1972).

Concentrada em análise de sobrevivência em vez deconfiança.

O modelo linear de taxa de falha proporcional assumeque

Cada variável x afeta a taxa de falha, oscilando para cima e para baixo.

A análise depende das pressuposições dos parâmetros

}exp{)();( 0 βTxthxth =

A análise depende das pressuposições dos parâmetros feitos em h0.

Se a taxa de falha de base acumulada , então:

onde y=logG(t).

ηαθ )]([);(0 tGtH =

)}exp()(exp{);( 0 βTxtHxtS −=

)}exp()]([exp{ βα η TxtG−=

)}logexp(exp{ βαη Txy ++−=


Segue que z = ηy + log α + xT β ~ Gumbele

y = -η-1 log α – xT(-η-1 β) + η-1z

sendo que esta é a fórmula do modelo log-linear deregressão da Weibull.

Pode-se usar outras fórmulas parah0 e a que possuir omelhor ajuste, pode então, ser utilizada.

5.2 Análise do Modelo Semiparamétrico deTaxa de Falha Proporcional

Deagoraemdiante,h0 nãoseráespecificado.Deagoraemdiante,h0 nãoseráespecificado.

Assume-sen itens,r tempos de falha distintos eRi é oconjunto de risco, ou seja, o grupo de itens que temchances de falhar um pouco antes det(i).

β é estimado pela função de verossimilhança:

∏∑=∈

=r

iRl

Tl

Ti

i

x

xL

1 )(

)(

)exp(

)exp()(

ββ

β


Há muitas justificativas para a verossimilhança acima:

1) Cox (1972) originalmente considerou análises condicionaisdos tempos de falha . A probabilidade condicional do item (i)falhar dado que conhecemos as falhas anteriores é:

A verossimilhança é formada pelo produto das probabilidadescondicionais.

(Basicamente, usa-se )

∑∑∈∈

=

ii Rl

Tl

Ti

Rlll

ii

x

x

xth

xth

)exp(

)exp(

);(

);(

)(

)(

)()(

)()(

ββ

∏=

−=r

iiir AAAAA

1111 ),...,Pr(),...,Pr(

2) Também, é obtida por Kalbfleisch & Prentice (1973) como averossimilhança

marginal deβ baseada nos postos das observações.

SeL(β ) é tratado como a verossimilhança, então a teoria damáximo verossimilhança pode ser aplicada:

i) O estimador de máximo verossimilhança é obtidomaximizandoL(β );

ii) Sua variância assintótica é obtida pela 2ª derivada del(β).

Computacionalmente:

Primeiro, calcula-se ;)exp( βλ Tl x=


Então, acumula-se o decréscimo dos tempos de falha pra obter , i = 1,...r.

L(β ) é formado pelo produto dos termos λ(i)/ci.

O log de verossimilhança e suas derivadas são:

Onde

∑∈

=iRl

lic λ

∑=

−=r

ii

Ti cxl

1)( )log()( ββ

∑=

−=r

iii vx

d

dl

1)( )(

)(ββ

∑=

−−=r

i

Tiii vvA

d

ld

12

2

)()(

ββ

Onde

e

A maximização numérica é usada para obter .

Há pacotes computacionais que podem nos auxiliar, porexemplo, o R.

Observações Empatadas:

Dados empatados complicam os cálculos da verossimilhança.

∑∈

=iRl

lli

i xc

v λ1∑∈

=iRl

Tlll

ii xx

cA λ1


Uma aproximação razoável neste caso é:

Ondedi é o número de tempos que falham emt(i) e s(i) éa soma das covariâncias que dependem destes tempos.

A aproximação é boa sedi /n(i) for pequeno.

As covariâncias que dependem dos tempos sãosubstituídas dex(l) porx(l) (t(i)).

5.3 Estimação de Sobrevivência e Função de Taxa de Falha

∏ ∑=∈

=r

iRl

dTl

Ti

i

ix

sL

1 )( )]exp([)exp(

)(β

ββ

Taxa de Falha

A taxa de sobrevivência de base é obtida por um método não paramétrico junto com o estimador do produto limite, como:

onde

Também, o estimador de H0 é dado por .

Os métodos do capítulo 3 podem ser aplicados aqui para verificar as formas paramétricas para S0 .

∏∑<

∈

−ttj

Rll

j

j

i

)(:^

^

)1(λ

λ )exp(^^

βλ ll x=

)(log)(^

0

^

0 tStH −=


5.4 Métodos de Verificação

As suposições da taxa de falha proporcional podem serverificadas pelos dados e modelos de ajuste:

Depois do ajuste, o gráfico do log vs. t deve ter paralelas verticais.

O resíduo pode ser definido por um item não censurado como:

)exp()();( 0 βTjj xthxth =

)exp()(^^

0 βTxtHe=

e por um item censurado, e deve ser somado por 1.

Se os resíduos se comportam como amostras de umadistribuição exponencial de média 1, o ajuste é adequado.

Deve-se ter cautela, pois os resíduos podem estaradequados mesmo quando a taxa de falha proporcionalnão estiver.

Também, o gráfico dos resíduos vs. as variáveis podenão se encaixar no modelo.


5.5 Exemplos Numéricos

Tabela 4.3 contem dados de falha sob stressFigura 4.7 indica a influencia da taxa de stress na falha

h(t;x) = h0(t)em(x) para alguma função m de x

caso m(x)= -β log x ,temos queβ^ =0,2069 com erro padrão assintótico =0,0457 => forte efeito da taxa de stress

Figura 5.1 mostra um gráfico de log^H0(t) vs log tindicandoWeibullindicandoWeibullA checagem da hipótese de modelo de taxas

proporcionais é feita segundo estratificação sugerida em§5.4

Isto leva ao gráfico na Figura 5.2 (paralelismo vertical=OK)

Checagem da linearidade em log x é feita através de um estimador não paramétrico dos efeitos de x

onde zj(x) é um indicador de que x está no nível j.Este modelo não é estimável (uma constante pode sempre ser adicionada ao αj)

( ) ( )∑=

=5

1ijj xzxm α


Mas se torna estimável ao se impor uma restrição (Σαj=0 ouαk =0, para algum k)

O valor maximizado de l é –177,00 em oposição a-178,19 sob o modelo linear.W = 2 (-177,00+178,19) = 2,38 < 7,81 =χ2

0,05(4-1) =>modelo linear aceito

Tabela 4.4 contem dados de tempo de falha de fibrasFigura 4.8 contem PH plots para diferentes níveis de

stress.Dando indicação de que o efeito dos lotes nos tempos

de falha do modelo proposto é

( ) ( )

+= ∑

7

logexp; zSthxth αγ

onde S é o nível de stress e zj é o indicador do lote j , j =1,..,7

( ) ( )

+= ∑=1

0 logexp;i

jj zSthxth αγ

Tabela 4.3 Teste de força em bilhas de liga de Tungstênio

0,1 1676 2213 2283 2297 2320 2412 2491 2527 2599 2693 2804 28611 1895 1908 2178 2299 2381 2422 2441 2458 2476 2528 2560 297010 2271 2357 2458 2536 2705 2783 2790 2827 2837 2875 2887 2899100 1997 2068 2076 2325 2384 2752 2799 2845 2899 2922 3098 31621000 2540 2544 2606 2690 2863 3007 3024 3068 3126 3156 3176 3685

Coluna 1: Taxa de Stress em MNm-2seg-1 ; Colunas 2-13: Tempo de falha


Figura 4.7Plot da Weibull para dados de tungestênio, com linhas ajustadas pela Weibull

Figura 5.1

*

*

**

**

**

**

*

7 . 4 7 . 5 7 . 6 7 . 7 7 . 8 7 . 9 8 . 0

-3-2

-10

1

L o g d o t e m p o d e f a l h a s o b a s t a x a s d e s t r e s s 0 . 1 ( * ) , 1 ( + ) , 1 0 ( - ) , 1 0 0 ( x ) , 1 0 0 0 ( o )

Lo

g(lo

g(p

rob

ab

ilid

ad

e d

e s

ob

revi

vênc

ia))

+

-

-

--

---

--

--

x

x

xx

xx

xx

xx

o

o

oo

ooo

o

*

****

* * * * **

7.4 7.6 7.8 8.0 8.2

-4-3

-2-1

01

2

Log do tempo de falha sob as taxas de stress 0.1(*),1(+),10(-),100(x),1000(o)

log

(Fun

ção

de

taxa

de

falh

a a

cum

ula

da

)

+

-

--

--

----

--

x

xx

xx

x xx xxx

o

oo

oo

ooooo

o


Figura 5.2

xx x xx

xx xxxxxxxxxxx

xxxxxxxxxxxx xxxxxxxxxx

xxxxxxxxxxxxxx

xx

7.6 7.8 8.0 8.2

-4-3

-2-1

01

23

log

(Fun

ção

de

taxa

de

falh

a a

cum

ula

da

)

Procedimentos iterativos funcionarão melhor casoiniciem com bons valores iniciais.

Estes podem ser obtidos através da MV com Weibullou por mínimos quadrados.γ^ =32,586 com erro padrão assintótico igual a 3,097 =>forte efeito da taxa de stress.

O teste para efeito do lote é baseado na Razão de MVentre M0: não há efeito do lote e M1: existe efeito do lote

A maximização das log verossimilhanças são –320,119e –261,775 que nos trazemW =2(-261,775+320119)=116,69 >>χ2

α(7)=> o efeito dolote é significante.

log(tempo de falha)


A hipótese de modelo de taxas proporcionais é testadapela estratificação dos dados através dos lotes.

A figura 5.3 mostra o gráfico do log^Hj aproximado porPH

A figura 5.4 mostra o gráfico do log^H0 vs log t , ondelinearidade indica Weibull

Também é possível investigar a interação entre o stresse os lotes.

Isto é feito através da inclusão de 7 covariáveis zj logS no modelo

A maximização da log-verossimilhança é –259,241 comW=5,07 (não significante)

Ajustando uma linha reta na figura 5.4 teremoslog^H0(t)= -11,19 + 1,533logt

EH^j(T)=15,296T1,533S32,586eα^j , j = 1,...,8 (α8 = 0)

Exemplo : o tempo médio de falha previsto para S= 23,4e j = 4 é T = 7,4 anos.


Tabela 4.4 Dados de Falha de fibra de Kevlar

Stress Lote Tempo Stress Lote Tempo Stress Lote Tempo Stress Lote Tempo29,7 2 2,2 29,7 5 243,9 27,6 2 694,1 25,5 1 11.487,329,7 7 4,0 29,7 4 254,1 27,6 4 876,7 25,5 5 11.727,129,7 7 4,0 29,7 1 444,4 27,6 1 930,4 25,5 4 13.501,329,7 7 4,6 29,7 8 590,4 27,6 6 1.254,9 25,5 1 14.032,029,7 7 6,1 29,7 8 638,2 27,6 4 1.275,6 25,5 4 29.808,029,7 6 6,7 29,7 1 755,2 27,6 4 1.536,8 25,5 1 31.008,029,7 7 7,9 29,7 1 952,2 27,6 1 1.755,5 23,4 7 4.000,029,7 5 8,3 29,7 1 1.108,2 27,6 8 2.046,2 23,4 7 5.376,029,7 2 8,5 29,7 4 1.148,5 27,6 4 6.177,5 23,4 6 7.320,029,7 2 9,1 29,7 4 1.569,3 25,5 6 225,2 23,4 3 8.616,029,7 2 10,2 29,7 4 1.750,6 25,5 7 503,6 23,4 5 9.120,029,7 3 12,5 29,7 4 1.802,1 25,5 3 1.087,7 23,4 2 14.400,029,7 5 13,3 27,6 3 19,1 25,5 2 1.134,3 23,4 6 16.104,0

Stress em Mpa; Numero do lote; Tempo de falha em horas(* caso dado censurado)

29,7 5 13,3 27,6 3 19,1 25,5 2 1.134,3 23,4 6 16.104,029,7 7 14,0 27,6 3 24,3 25,5 2 1.824,3 23,4 5 20.231,029,7 3 14,6 27,6 3 69,8 25,5 2 1.920,1 23,4 6 20.233,029,7 6 15,0 27,6 2 71,2 25,5 2 2.383,0 23,4 5 35.880,029,7 3 18,7 27,6 3 136,0 25,5 3 2.442,5 23,4 1 41.000,029,7 2 22,1 27,6 2 199,1 25,5 8 2.974,6 23,4 1 41.000,029,7 7 45,9 27,6 2 403,7 25,5 2 3.708,9 23,4 1 41.000,029,7 2 55,4 27,6 2 432,2 25,5 8 4.908,9 23,4 1 41.000,029,7 7 61,2 27,6 1 453,4 25,5 2 5.556,0 23,4 4 41.000,029,7 5 87,5 27,6 2 514,1 25,5 6 6.271,1 23,4 4 41.000,029,7 8 98,2 27,6 6 514,2 25,5 8 7.332,0 23,4 4 41.000,029,7 3 101,0 27,6 6 541,6 25,5 8 7.918,7 23,4 4 41.000,029,7 2 111,4 27,6 2 544,9 25,5 6 7.996,0 23,4 8 41.000,029,7 6 144,0 27,6 8 554,2 25,5 8 9.240,3 23,4 8 41.000,029,7 2 158,0 27,6 1 664,5 25,5 8 9.973,0 23,4 8 41.000,0


Figura 4.8 plot da Weibull para dados de fibra de Kevlar sob pressão , com linhas ajustadas pela Weibull

Figura 5.3Log da função de taxa de falha acumulada

+

++ ++

+++++

+++++++ +++ ++++

++ ++ + ++++++ ++

2 4 6 8 1 0 1 2

-5-4

-3-2

-10

1

L o g d o te m p o d e fa lha no s níve is d e s tre s s 2 9 .7 (+ ),2 7 .6 (-),2 5 .5 (x),2 3 .4 (o )

Lo

g(-

log

(pro

ba

bili

da

de

de

so

bre

vivê

ncia

))

-

--- - - --

------

-- -- ----

-

x

xxx xx

xxx x xxx

xxxxxx

xxx x

o

oooo

oooo

o

Figura 5.3Log da função de taxa de falha acumulada para dados da fibra de Kevlar: modelo estratificado . Lotes 1-8

1

111 11

1 1 1 1 1

0 2 4 6 8 1 0 1 2

-6-4

-20

24

lo g ( te m p o d e fa lh a )

log

(Fun

ção

de

taxa

de

falh

a a

cum

ula

da

)

2222 2 2 2 2 2 2 22 22 2 2 22

2 2 2

3

333

3 3 3 33

3

4

444 44

444 4 4

5

55

55

55

6

66

666

6 6666

77

77 7 7 7

77

8

888

8 8 8 8888


Figura 5.4 Log da Log da função de taxa de falhaacumulada para dados da fibra de Kevlar

o

oo o

o oooo o ooooo oo o o o o o oo o oo o oo o o o oo oooo oooooo

o ooo o ooo oooo oo oo oooooo oo o oo o oo oo oooo oooo oo ooo o oo oo o

0 2 4 6 8 1 0 1 2

-10

-50

5

lo g (te m p o d e fa lha )

log

(Fun

ção

de

taxa

de

falh

a a

cum

ula

da

)

Em geral:

i = 0, 1, ..., n-2


6. MODELOS BAYESIANOS DINÂMICOS

6.1 Introdução

Os modelos PH (Proportional Harzard Model) são razoáveisem alguns casos, porem muitas as vezes modelos maisgenéricos são requeridos.

Os enfoques utilizados são :

1 - Parametrização das mudanças entre as faixas de estimação. Exemplo :

2 - Estimação independente por partes em cada

kte−= 0ββ2 - Estimação independente por partes em cada

intervalo.

O primeiro enfoque exige uma suposição forte em tornodos dados enquanto o segundo não exige nada ,porem asinformações contidas nos intervalos anteriores sãodesconsideradas.

Modelos dinâmicos permite mapearmos as mudanças entreas diferentes faixas sem a necessidade de fortes suposiçõesepermitindo que os dados orientem a direção do modelo.


Modelos dinâmicos provem uma abordagem para osmodelos de PHM.

6.2 Elementos do Modelo

Distribuição Exponencial por Partes

A Distribuição Exponencial por Partes é usada paratempo de falha baseada na função de densidade Harzard :

, onde o primeiro componente (baseline) e o primeiro componente .

iTx

i exth βλ ==);( ..2,1, =∈ iIt i

0,ii BB => 1=tx

Fatorização Temporal da verossimilhança

Os dados podem ser divididos ao longo do mesmointervalo Ii . Cada intervalo prove a verossimilhançabaseado na condicional das informações obtidos dointervalo anterior.


Esta é a verossimilhança do intervalo Ii condicionado em Di-1,F(Ii) é o conjunto de itens que falharam no intervalo Ii ,R(Ii) é o conjunto de itens que estavam em risco antes Ii mais não falharam em Ii eL é o conjunto de parâmetros

6.3 Desenvolvendo a Análise

∏∏

∏

=−

=−

=

∧∧=

=

N

Iiiq

N

Iiwi

N

Ii

DtSDtpL

onde

LLTFL

11

11

1

),|(),|(

:

6.3 Desenvolvendo a Análise

Evolução Paramétrica

Como descrito em 4.6, parâmetros sucessivos podem ser calculados via

WbW

e

Itamanhobi

WwVar

wE

onde

wGi

ii

i

ii

i

iiii

=

==

=

+−=

)(

)(

0)(

,

)( 1βββ


O valor de W indica o potencial de mudança dos b´spor unidade de tempo.

W ⇒ 0 conduz ao um modelo estático.W-1 ⇒ 0 grande potencial de mudança conduzindo amodelos independentes.

Assume-se que : (especificação parcial)

A priori para o intervalo Ii é onde

],[~| 11 −− iiii CmDβ

],[~| 1 iiii PD αβ −

iii mbGi= −)( 1α

Novamente, o valor de Wi (ou W) podem ser definidos como fatores de desconto.

Verossimilhança e Atualização

A verossimilhança Li pode ser escrita como :

iiT

iiii WbGCbGP += − )()( 1

∏∏∏=∈

−−

∈

−− == −−i

i

iiqq

i

i

iww

i

r

jij

IRq

tt

IFw

ttwi LeeLi

1)(

))((

)(

))(()( 1)(

1)( λλλ


Onde

Onde é o valor constante da hazard para o item s em Ii , ri é o numero do item em Iidij é o indicador de falha para o item em j ,tij é o tempo de sobrevivência do indivíduo j em Ii .

Uma analise baysiana completa envolve a especificaçãodaprior para atéa integraçãonuméricade cada

))(()(

1)(

))((

)(

))(()(

1)(

1)(

1)(

][ −

−−

−−

=∈

−−

∈

−−

=

== ∏∏∏

iijj

iij

i

i

iiqq

i

i

iww

i

ttjiij

r

jij

IRq

tt

IFw

ttwi

eL

LeeLi

λδ

λλ

λ

λ

βλ)()( T

sxsi e=

daprior para atéa integraçãonuméricade cadaintervalo Ii , alternativamente, especificações parciaispodem ser utilizadas.A atualização é efetuada em ciclos através dos intervalos

e dentro de cada intervalo através de seus itens. O calculoenvolvendo esta analise esta especificado abaixo:

oD|1β

1,1~

| −− jii Dβ

ijj

i ZRG βλ~

)(log:.. ≅

1,1)( | −− ji

ji Dλ

ijl

jij

i D ,1)( | −λ

jii D ,1~

| −β1+⇒ jj

jij

ii D ,1)(

~,| −λβ

LINEAR BAYES

MEDIAiλ

SAIDARodrigues Fazenda 122

Passo 1 : Dado obtemos :

Passo 2 : Aproximar por uma Gammautilizando os dados do Passo 1. Parâmetros

e ;

Passo3 : Atualizar, usandoa informaçãodo Lij para

],[~| 1;1 ijijjii CD αβ −−

;

,

],[~|log 1,1)(

jijT

jij

ijT

jij

ijijjij

i

xCxq

axf

onde

qfD

=

=

−−λ

1,1)( | −− ji

ji Dλ

1)( −= ijij qα )(1)( ijfij eqij

−−=γ

Passo3 : Atualizar, usandoa informaçãodo Lij para;

Passo 4 : Use os métodos lineares de bayes e a relaçãoentre e para obter

Passo 5 : Defina ai,j+1 = mij e Pi,j+1=Cij e repita ospasso anteriores;

Passo 6 : Quando j=rij , todos os itens Li tem sidoprocessado. Evolui para o próximo intervalo.

],[~| 1,1)(

−− −++ iijijijijjij

i ttGD γδαλ

)(log jiλ iβ

];,[~| ,1~

ijijjii CmD −β


O procedimento acima tem como resultado aespecificação parcial da distribuição :

Este procedimento continua ate In , o ultimo intervalocom informações.

SmoothingÉ mais conveniente trabalhar com a distribuição de

Ela gera resultados mais precisos e é utilizada napredição.Existem distribuições Smoothing obtidas via

,...;2,1,|~

=iDiiβ

.,..,1,|~

NiDNi =β

predição.Existem distribuições Smoothing obtidas viaalgoritmos recursivos.

Com valores iniciais

iiTii

Niiii

Tiii

Ni

iNiii

Tiii

Ni

Ni

Nini

CbGPCPPbGCCC

amPbGCmm

então

CmD

)()()(

)()(

],[~|

11111

111

11111

~

+−+++

−++

++−++

−−=

−+=

β

NNN mm = N

NN CC =


6.4 Predição de Ocorrências Futuras

Dois casos devem ser considerados :

1) Predição dos itens que sobraram e ainda não falharam2) Predição de novos itens.Para o caso 2, a distribuição de sobrevivência predita éS(W\DN) onde ,

i

i

jNjjNiN

Iw

com

DtWtSDtWwSDWS

∈

>>= ∏−

=−−

1

111 ),|(),|()|(

O fato acima pode ser obtido da seguinte forma :

Obtido de

Aproximados pela distribuição de .

i

∫

∫∞

−−−

∞−−−

>=

>>=>

−0 1

))((

0 111

),|(

),|(),,|(),|(

1kNkk

tw

kNkkNkkNi

dDtWpe

dDtWpDtWwSDtWwS

kk λλ

λλλ

λ

),|( 1 Nkk DtWp −>λ ),|( 1 Nkk DtWp −>β

Nk D|β


Se então,

Cálculos parecidos são utilizados para a distribuiçãopreditora a posteriori .

Devemos verificar se a especificação da priori e razoável.

6.5 Ilustrações

),(~| kkNk GD γαλ

i

i

j i

i

i

iN

Iw

btwDwS ii

∈

+−+= ∏−

=

−−1

1

)1()1()|( αα

γγ

)|( 0DwS

6.5 Ilustrações

Exemplo 1 : Dados de 90 pacientes com câncer gástricodivididos em dois grupos de tratamento : Quimioterapia eRadioterapia / somente quimioterapia.

PHM ajustado com 2 covariáveis : x1 - Indicador de tratamento e x2=x1*t. Isto é equivalente a assumir um modelo de regressão com os coeficientes .

Os estimadores MLE são :

t211 Θ+Θ=β

2711.1^

1 =Θ 0794.0^

2 =Θ


Os estimadores MLE são :

Analises dinâmicas foram feitas com :

1 - Priori Vaga

2 - dados obtidos pelos tempos observados de falha.

3 - evolução do erro da variância , dado por :

2711.1^

1 =Θ 0794.0^

2 =Θ

];10,0[~| 23

01 IDβ

=0016.00

00biWi

Para a variação foi utilizado um crescimento linearcom em 4.6 .

Figura 1 mostra as estimativas do coeficiente deregressão.

1β


Figura 2 mostra a curva de sobrevivência predita.Figura 2 mostra a curva de sobrevivência predita.


Figura 3 mostra estimativas para as taxas de mudança de .1β

Exemplo 2 : Dados de 300 desempregados que foramacompanhados durante um ano na Inglaterra. Esta é umaamostra de um grande estudo feito por Germerman eWest (1987). As covariáveis selecionadas são x1 ~ idadee x2 ~ faixa de renda.Analises dinâmicas foram feitas com :

1 - Priori

2 - Intervalos Semanais;

3 - Evolução do erro da variância Wi=7*diag(10-2,10-5,10-3).

)];)025.0(,)0125.0(,1(,)25.0,025.0,1[(~| 2201 diagD T−−−β


Figura 4 mostra preposteriori par 4 hipotéticosindivíduos.

Figura 5 mostra estimativasparamétricas(prioriestavaFigura 5 mostra estimativasparamétricas(prioriestavamuito forte).


Figura 6 mostra predição para os mesmos 4 indivíduos.


Documents

Análise de Sobrevivência€¦ · Análise de Sobrevivência é o estudo de indivíduos (itens observados) onde um evento bem definido (falha) ocorre depois de algum tempo (tempo