disciplinas:ce089-2014-02

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » disciplinas » ce089-2014-02

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- disciplinas:ce089-2014-02 [2014/08/18 19:55]
walmes
+++ disciplinas:ce089-2014-02 [2014/12/15 14:16] (atual)
walmes
@@ Linha 22: / Linha 22: @@
 ----
+/*
 ==== Histórico das Aulas do Curso ====
@@ Linha 37: / Linha 37: @@
   - 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}}.
 ----
@@ Linha 46: / Linha 207: @@
 ==== Avaliações ====
-  * Trabalho 1: {{Geração de números aleatórios|http://www.leg.ufpr.br/~walmes/ensino/ce089-2014-02/ce089-2014-02-trab01.html}};
+  * 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}}

Diferenças

Navegação

Busca

Ferramentas

QR Code