CE-083: Estatística Computacional II

"The purpose of computing is insight, not numbers." (Richard Hamming)

• Professor: Walmes Marques Zeviani, (LEG: Laboratório de Estatística e Geoinformação)
• Curso: Estatística.
• Período: 2014/2.
• Local: LABEST LAB C.
• Horário: Segunda 20h45-22h00 e Quarta 19:00-20:30h.
• Atendimento: Segunda, 19:00-20:30h.

##-----------------------------------------------------------------------------
 
require(bbmle)
 
##-----------------------------------------------------------------------------
 
y <- rnorm(100, mean=0, sd=1)
crt <- 1.5
r <- y<crt
y[!r] <- crt
r <- as.integer(r)
plot(ecdf(y))
table(r)
 
ll <- function(theta, y, r){
    ys <- split(y, r)
    ll1 <- dnorm(ys[["1"]],
                 mean=theta[1],
                 sd=exp(theta[2]),
                 log=TRUE)
    ll2 <- pnorm(ys[["0"]],
                 mean=theta[1],
                 sd=exp(theta[2]),
                 lower.tail=FALSE,
                 log=TRUE)
    ll12 <- sum(ll1)+sum(ll2)
    ## Tem que retornar o negativo da soma. A mle2() minimiza.
    return(-ll12)
}
 
##-----------------------------------------------------------------------------
## Estimação com a bbmle.
 
parnames(ll) <- c("mu", "lsigma")
m0 <- mle2(minuslogl=ll, start=list(mu=0, lsigma=log(1)),
           vecpar=TRUE, data=list(y=y, r=r), method="BFGS")
 
coef(m0)
summary(m0)
ci <- confint(m0, method="quad")
ci
 
htheta <- coef(m0)
logL <- c(logLik(m0))
 
##-----------------------------------------------------------------------------
 
llv <- Vectorize(
    FUN=function(th1, th2, y, r){
        -ll(c(th1, th2), y=y, r=r)
    },
    vectorize.args=c("th1", "th2"))
 
th1 <- seq(-0.4, 0.4, l=30)
th2 <- seq(-0.1, 0.4, l=30)
 
lla <- outer(th1, th2, llv, y=y, r=r)
 
contour(th1, th2, lla,
        xlab=expression(mu),
        ylab=expression(log(sigma)))
abline(v=htheta[1], h=htheta[2], lty=2)
contour(th1, th2, lla, add=TRUE,
        levels=(logL-0.5*qchisq(0.95, df=length(htheta))),
        col=2)
abline(v=ci[1,], h=ci[2,], lty=3, col=2)
 
plot(profile(m0))
 
 
##--------------------------------------------
## Implementação conforme sugestão do Wikipedia.
 
require(rootSolve)
 
n <- rbinom(1, size=200, p=0.8)
y <- c(rpois(n, lambda=exp(2)), rep(0, 200-n))
length(y)
 
barx <- mean(y)
p <- sum(y==0)/length(y)
 
f <- function(lambda, L){
    L$barx*(1-exp(-lambda))-lambda*(1-L$p)
}
 
L <- list(barx=barx, p=p)
gradient(f, x=2, L=L)
 
curve(f(x, L=L), 0, 15); abline(h=0)
 
## Newton-Raphson.
maxiter <- 50; i <- 1          ## Número máximo de iterações e contador.
tol <- 1e-5; error <- 100*tol  ## Tolerância e erro inicial.
theta <- matrix(NA, nrow=1, ncol=1)
theta[1,] <- barx
while(i <= maxiter & error>tol){
    theta <- rbind(theta, rep(NA,1))
    G <- f(theta[i,], L)
    H <- gradient(f=f, x=theta[i,], L=L)
    theta[i+1,] <- theta[i,]-solve(H)%*%G
    error <- sum(abs((theta[i+1,]-theta[i,])/theta[i+1,]))
    i <- i+1
    ## print(c(theta[i,], G))
    print(cbind(H, G, theta[i,]))
}
 
lam <- theta[i,]        ## Estimativa do lambda.
pii <- 1-barx/theta[i,] ## Estimativa do pi.
 
##-----------------------------------------------------------------------------
## Vendo os contornos da verossimilhança.
 
llmax <- ll(th=c(log(pii/(1-pii)), log(lam)), y=y)
 
th1 <- seq(-7, 5, l=50)
th2 <- seq(-1, 4, l=50)
lla <- outer(th1, th2, llv, y=y)
 
contour(th1, th2, lla,
        ## levels=seq(from=llmax-30, to=llmax, by=10),
        nlevels=20,
        xlab="theta1: logit(p)",
        ylab="theta2: log(lambda)")
abline(v=log(pii/(1-pii)), h=log(lam), col=2)

1. Instruções de como enviar os trabalhos/atividades: datafilehost.

• Trabalho 1: Geração de números aleatórios;