## Modelo Binomial Negativo ## Gerando a amostra set.seed(123) y <- rnbinom(100, size = 100, p = 0.8) ## Vendo os dados hist(y) ## Passo 1 - Função de log-verossimilhança ll.neg.binomial <- function(par, y) { phi <- exp(par[1]) p <- par[2] n <- length(y) ll <- n * phi * log(p) + sum(y) * log(1 - p) + sum(log(gamma(y + phi))) - n * log(gamma(phi)) - sum(log(factorial(y))) #print(c(phi,p)) return(ll) } ## avaliando a log-verossimilhança no ponto ll.neg.binomial(par=c(log(100),0.8),y=y) ## Precisamos maximizar a log-verossimilhança ## Maximizando pela optim estimativa.bin <- optim(par=c(log(50),0.5), fn = ll.neg.binomial, y = y, method="L-BFGS-B", lower=c(1e-32,1e-32), upper=c(Inf, 1), control=list(fnscale=-1), hessian=TRUE) estimativa.bin ## Construindo o intervalo de confiança de Wald V.esti <- solve(-estimativa.bin$hessian) std.error <- sqrt(diag(V.esti)) ic.phi <- estimativa.bin$par[1] + c(-1,1)* qnorm(0.975)*std.error[1] ic.phi ic.p <- estimativa.bin$par[2] + c(-1,1)* qnorm(0.975)*std.error[2] ic.p ## Perfil de Verossimilhança p.gride <- seq(0.5,0.85,l=100) perfil.p <- function(phi,p, y) { n <- length(y) ll <- n * phi * log(p) + sum(y) * log(1 - p) + sum(log(gamma(y + phi))) - n * log(gamma(phi)) - sum(log(factorial(y))) print(c(phi,p)) return(ll) } profile.p <- c() for(i in 1:100){ profile.p[i] <- optim(par = 50, fn = perfil.p, p=p.gride[i], y=y, method="L-BFGS-B", lower=c(1e-32), upper=c(130), control=list(fnscale=-1))$value print(i) } dev <- 2*(estimativa.bin$value - profile.p) plot(dev ~ p.gride,type="l",xlim=c(0,1)) abline(h=3.84) abline(v=ic.p) ########################################################################################## ## Otmizando usando o gradiente analitico ################################################ ########################################################################################## ## programar a função escore escore <- function(par, y) { phi <- par[1] p <- par[2] n <- length(y) U.phi <- n * log(p) + sum(digamma(y + phi)) - n * digamma(phi) U.p <- (n * phi)/p - sum(y)/(1 - p) return(c(U.phi, U.p)) } ## Usando optim com escore analitico escore(par=estimativa.bin$par,y=y) system.time(est.bin <- optim(par=c(25,0.5), fn = ll.neg.binomial, y = y, method="L-BFGS-B", lower=c(1e-32,1e-32), upper=c(Inf, 1), control=list(fnscale=-1), hessian=TRUE)) system.time(est.bin.esc <- optim(par=c(25,0.5), fn = ll.neg.binomial, y = y, gr= escore, method="L-BFGS-B", lower=c(1e-32,1e-32), upper=c(Inf, 1), control=list(fnscale=-1), hessian=TRUE)) est.bin.esc$par est.bin$par est.bin.esc$hessian est.bin$hessian ### Usar o algoritmo de NewtonRaphson que exige a segunda derivada Hessiana <- function(par, y) { phi = par[1] p = par[2] n <- length(y) II <- matrix(c(sum(trigamma(y + phi)) - n * trigamma(phi), n/p, n/p, -(n * phi)/p^2 - sum(y)/(1 - p)^2), 2, 2) return(II) } Hessiana(par=est.bin.esc$par,y=y) est.bin.esc$hessian ## Usando o Newton Raphson NewtonRaphson <- function(initial, escore, hessiano, tol = 1e-04, max.iter, n.dim, ...) { solucao <- matrix(NA, max.iter, n.dim) solucao[1, ] <- initial for (i in 2:max.iter) { solucao[i, ] <- initial - solve(hessiano(initial, ...), escore(initial, ...)) initial <- solucao[i, ] tolera <- abs(solucao[i, ] - solucao[i - 1, ]) # print(initial) if (all(tolera < tol) == TRUE) break } return(initial) } ## Aplicando na função Binomial Negativa system.time(est.nr <- NewtonRaphson(initial = c(50, 0.5), escore = escore, hessiano = Hessiana, max.iter = 100, tol = 1e-10, n.dim = 2, y = y)) est.nr est.bin.esc$par est.bin$par ## Encontrar uma região de confiança baseada na Deviance p <- seq(0.5, 0.9, l=50) phi <- log(seq(10,130,l=50)) gride <- expand.grid(phi,p) plot(gride) dim(gride) ## Avaliando a função a log-verossimilhança ll <- c() for(i in 1:2500){ ll[i] <- ll.neg.binomial(par= as.numeric(gride[i,]), y=y) } ll dev <- 2*(est.bin.esc$value - unlist(ll)) contour(x=p, y = phi, z = matrix(dev,50,50), levels= qchisq(seq(0.1,0.99,l=5),df=2)) ## Distribuição empirica est.nr <- matrix(NA, ncol=2, nrow=500) for(i in 1:500){ y <- rnbinom(100, size=10, p =0.8) temp <- try(optim(par=c(log(10),0.5), fn = ll.neg.binomial, y = y, control=list(fnscale=-1), hessian=TRUE)$par) if(class(temp) != "try-error"){ est.nr[i,] <- temp} } par(mfrow=c(1,2)) hist(est.nr[,1]) hist(est.nr[,2])