#
#Script para estimativa dos parâmetros a (theta ou alfa) e b (beta) do 
# PLP, tempo ótimo entre preventivas, e intervalos de confiança
#
 estimatotpl<-function(tij,b,II,Lambda,sign,Cmc,Cmp)
 {
#Estimador MVS (lei de potencias) para i sistemas com j falhas#
#                                                            #
############  tij é a matrix dos tempos de falhas  ############
#                                                            #
	Delta <- 0.1
	itera <- 0
	while(itera < II & abs(Delta) > Lambda) {
		NN <- length(tij[, 1])	# comprimento da matriz tij #
		N <- 0	# Numero de ocorrencias (exeto censuras) #
		STB <- 0	# Soma dos t's elevado a b #
		STBLT <- 0	# Soma dos t's elevado a b vezes o log do t. #
		SLTij <- 0	# Soma do log dos t's #
		for(i in 1:(NN - 1)) {
			if(tij[i, 3] == 1) {
				N <- N + 1
				SLTij <- SLTij + log(tij[i, 2])
			}
			if(tij[i, 1] != tij[(i + 1), 1]) {
				STB <- STB + (tij[i, 2])^b
				STBLT <- STBLT + ((tij[i, 2])^b) * log(tij[i, 2])
			}
		}
		STB <- STB + (tij[NN, 2])^b
		STBLT <- STBLT + ((tij[NN, 2])^b) * log(tij[NN, 2])
		Delta <- abs(b - N/((N/STB) * STBLT - SLTij))
		b <- N/((N/STB) * STBLT - SLTij)
		a <- (STB/N)^(1/b)
		itera <- itera + 1
	}
#                                                                                                                   #
#  CONSTRUÇÃO DE INTERVALO DE CONFIANÇA PARA BETA  #
#                                                  #
	STB <- 0
	STBLT <- 0
	STBLTa <- 0
	STBLTLTa <- 0
	for(i in 1:(NN - 1)) {
		if(tij[i, 1] != tij[(i + 1), 1]) {
			STB <- STB + (tij[i, 2])^b
			STBLT <- STBLT + ((tij[i, 2])^b) * log(tij[i, 2])
			STBLTa <- STBLTa + ((tij[i, 2])^b) * log(tij[i, 2]/a)
			STBLTLTa <- STBLTLTa + ((tij[i, 2])^b) * log(tij[i, 2]) * log(tij[i, 2]/a)
		}
	}
	STB <- STB + (tij[NN, 2])^b
	STBLT <- STBLT + ((tij[NN, 2])^b) * log(tij[NN, 2])
	STBLTa <- STBLTa + ((tij[NN, 2])^b) * log(tij[NN, 2]/a)
	STBLTLTa <- STBLTLTa + ((tij[NN, 2])^b) * log(tij[NN, 2]) * log(tij[NN, 2]/a)
	D2a2 <- (b * N)/(a^2) - (b * (b + 1) * STB)/(a^(b + 2))
	D2b2 <-  - N/(b^2) + (log(a) * STBLTa)/(a^b) - STBLTLTa/(a^b)
	D2a1b1 <-  - N/a + (STB - b * log(a) * STB + b * STBLT)/(a^(b + 1))
	MATRIZ <- matrix(nrow = 2, ncol = 2)
	MATRIZ[1, 1] <- D2a2
	MATRIZ[1, 2] <- D2a1b1
	MATRIZ[2, 1] <- D2a1b1
	MATRIZ[2, 2] <- D2b2
	INFOR <- ( - solve(MATRIZ))
	LS <- b + qnorm(1 - (100 - sign)/200) * sqrt(INFOR[2, 2])
	LI <- b - qnorm(1 - (100 - sign)/200) * sqrt(INFOR[2, 2])
	cat("Beta=", b, "  (", LI, ",", LS, ")", sign, "%", "\n")
	cat("Alfa=", a, "\n")
	cat("Var(Alfa) = ", INFOR[1, 1], "Var(Beta) = ", INFOR[2, 2], "Cov(Alfa, Beta) = ", INFOR[2, 1], "\n")
#
#  CALCULO DO TEMPO OTIMO  #
#                          #
	tot <- a * (Cmp/((b - 1) * Cmc))^(1/b)	#  INTERVALO DE CONFIANÇA PARA O TEMPO OTIMO  #
#                           #
	Dtotda <- ((Cmp/Cmc) * (1/(b - 1)))^(1/b)
	Dtotdb <- (a * (Cmp/((b - 1) * Cmc))^(1/b)) * (log(Cmp/((b - 1) * Cmc))/b^2 + 1/(b * (b - 1)))
	Vartot <- INFOR[1, 1] * (Dtotda^2) + INFOR[2, 2] * (Dtotdb^2) + 2 * INFOR[1, 2] * Dtotda * Dtotdb
	Li <- tot - qnorm(1 - (100 - sign)/200) * (Vartot^(1/2))
	Ls <- tot + qnorm(1 - (100 - sign)/200) * (Vartot^(1/2))
	cat("Tempo otimo entre preventivas  tot=", tot, "(", Li, ",", Ls, ")", sign, "%", "\n")
  }
#
#
 tij<-read.table("tij.txt",header=TRUE)
 b<-5
 II<-500
 Lambda<-1e-005
 sign<-95
 Cmc<-15
 Cmp<-1
 estimatotpl(tij,b,II,Lambda,sign,Cmc, Cmp)
#