====== CE-083: Estatística Computacional II ======

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

{{ http://www.leg.ufpr.br/~walmes/img/cupCoffeeRefill.jpg?400|}}

==== Detalhes da oferta da disciplina ====

  * Professor: [[http://www.leg.ufpr.br/~walmes|Walmes Marques Zeviani]], ([[http://www.leg.ufpr.br|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.

----

==== Scripts, notas e documentos ====

{{url>http://www.leg.ufpr.br/~walmes/ensino/ce089-2014-02/ 800px, 600px center}}

----

/*
==== Histórico das Aulas do Curso ====

{{ http://www.leg.ufpr.br/~walmes/img/tumblr_mvewnudGLK1r55ed0o2_400.gif?200|}}


  - 04/08:
    * Boas vindas;
    * Detalhes do Curso.
  - 06/08:
    * Geração de números aleatórios;
    * Método da congruência;
    * Método da transformação integral da probabilidade;
  - 18/08:
    * Como escrever documentos R+Markdown, por Vanessa Sehaber e Henrique Laureano.
  - 20/08:
    * Método de aceitação e rejeição.
  - 25/08:
    * Método Box-Muller para gerar números aleatórios da distribuição Gaussiana.
  - 27/08:
    * Métodos de geração de números aleatórios baseados em MCMC.
    * Amostrador independente (independence sampler).
    * Metropolis Random-Walk.
  - 01/09:
    * Teste de hipótese Monte Carlo.
  - 03/09:
    * Avaliação do nível de significância nominal de testes de hipótese.
  - 10/09:
    * Nível de cobertura de intervalos de confiança.
  - 15/09:
    * Avaliação via simulação do poder de testes paramétricos e não paramétricos na presença e ausência dos pressupostos.
  - 17/09:
    * Avaliação via simulação do poder de testes paramétricos e não paramétricos na presença e ausência dos pressupostos.
  - 22/09:
    * Nível de cobertura de intervalos de confiança - intervalos baseados na deviance a na sua aproximação (Wald).
  - 24/09:
    * Nível de cobertura de intervalos de confiança - efeito da parametrização.
  - 29/09:
    * Desenvolvimento do trabalho 03.
  - 01/10:
    * Desenvolvimento do trabalho 03.
  - 13/10:
    * Método Newton-Raphson.

----

*/

<code R>
##-----------------------------------------------------------------------------

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)

</code>

==== Links úteis ====

  - Instruções de como enviar os trabalhos/atividades: {{http://www.leg.ufpr.br/doku.php/disciplinas:datafilehost}}.

----

==== Avaliações ====

  * Trabalho 1: {{http://www.leg.ufpr.br/~walmes/ensino/ce089-2014-02/ce089-2014-02-trab01.html|Geração de números aleatórios}};

{{threads>pessoais:walmes:ce089-2014-02:discussion}}

~~DISCUSSION~~