Skip to content
Snippets Groups Projects
Commit 637fd5c9 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Remove a chamada de pacotes no meio.

parent fc435800
Branches
No related tags found
No related merge requests found
...@@ -81,14 +81,13 @@ gamma <- 0 ...@@ -81,14 +81,13 @@ gamma <- 0
fy <- dpg0(y = y, theta = theta, gamma = gamma) fy <- dpg0(y = y, theta = theta, gamma = gamma)
plot(fy ~ y, type = "h") plot(fy ~ y, type = "h")
lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2) lines(y + 0.3, dpois(y, lambda = theta), type = "h", col = 2,
xlab = "y", ylab = "f(y)")
``` ```
## Recursos interativos com o `rpanel` ## ## Recursos interativos com o `rpanel` ##
```{r, eval=FALSE} ```{r, eval=FALSE}
library(rpanel)
react <- function(panel){ react <- function(panel){
with(panel, with(panel,
{ {
...@@ -213,7 +212,7 @@ rp.checkbox(panel = panel, ...@@ -213,7 +212,7 @@ rp.checkbox(panel = panel,
rp.do(panel = panel, action = react) rp.do(panel = panel, action = react)
``` ```
## O espaço paramétrico ## ## O Espaço Paramétrico ##
```{r} ```{r}
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
...@@ -224,16 +223,18 @@ rp.do(panel = panel, action = react) ...@@ -224,16 +223,18 @@ rp.do(panel = panel, action = react)
# dpg1(y = 0:10, lambda = 1, alpha = -0.1) # dpg1(y = 0:10, lambda = 1, alpha = -0.1)
# undebug(dpg1) # undebug(dpg1)
library(latticeExtra)
y <- 0:200
fun <- Vectorize(vectorize.args = c("theta", "gamma"), fun <- Vectorize(vectorize.args = c("theta", "gamma"),
FUN = function(theta, gamma) { FUN = function(theta, gamma) {
sum(dpg0(y = y, theta = theta, gamma = gamma)) sum(dpg0(y = y, theta = theta, gamma = gamma))
}) })
grid <- list(theta = seq(0.5, 50, by = 0.5), grid <- list(theta = seq(1, 50, by = 1),
gamma = seq(-0.98, 0.98, by = 0.02)) gamma = seq(-0.5, 1, by = 0.05))
str(grid)
y <- 0:500
my <- max(y)
grid$sum <- with(grid, outer(theta, gamma, fun)) grid$sum <- with(grid, outer(theta, gamma, fun))
grid <- with(grid, grid <- with(grid,
cbind(expand.grid(theta = theta, gamma = gamma), cbind(expand.grid(theta = theta, gamma = gamma),
...@@ -242,16 +243,16 @@ grid <- with(grid, ...@@ -242,16 +243,16 @@ grid <- with(grid,
levelplot(sum ~ theta + gamma, levelplot(sum ~ theta + gamma,
data = subset(grid, round(sum, 3) == 1), data = subset(grid, round(sum, 3) == 1),
col.regions = gray.colors) + col.regions = gray.colors) +
layer(panel.abline(a = 0, b = -1/200)) layer(panel.abline(a = 0, b = -1/my)) +
layer(panel.abline(h = 0, lty = 2))
#-----------------------------------------------------------------------
fun <- Vectorize(vectorize.args = c("lambda", "alpha"), fun <- Vectorize(vectorize.args = c("lambda", "alpha"),
FUN = function(lambda, alpha) { FUN = function(lambda, alpha) {
sum(dpg1(y = y, lambda = lambda, alpha = alpha)) sum(dpg1(y = y, lambda = lambda, alpha = alpha))
}) })
# dpg1(y = 0:10, lambda = 5, alpha = -0)
# dpois(0:10, lambda = 5)
grid <- list(lambda = seq(0.2, 50, by = 0.2), grid <- list(lambda = seq(0.2, 50, by = 0.2),
alpha = seq(-0.98, 0.98, by = 0.02)) alpha = seq(-0.98, 0.98, by = 0.02))
grid$sum <- with(grid, outer(lambda, alpha, fun)) grid$sum <- with(grid, outer(lambda, alpha, fun))
...@@ -262,15 +263,16 @@ grid <- with(grid, ...@@ -262,15 +263,16 @@ grid <- with(grid,
levelplot(sum ~ lambda + alpha, levelplot(sum ~ lambda + alpha,
data = subset(grid, round(sum, 3) == 1), data = subset(grid, round(sum, 3) == 1),
col.regions = gray.colors) col.regions = gray.colors) +
layer(panel.abline(h = 0, lty = 2)) +
layer(panel.curve(-1/x))
# Já que lambda * alpha > -1, então alpha = -1/lambda dá a fronteira.
``` ```
## Verossimilhança e Estimação ## ## Verossimilhança e Estimação ##
```{r} ```{r}
library(bbmle)
library(corrplot)
# Função de log-Verossimilhança da Poisson Generalizada na # Função de log-Verossimilhança da Poisson Generalizada na
# parametrização de modelo de regressão. # parametrização de modelo de regressão.
llpg <- function(theta, y, X, offset = NULL) { llpg <- function(theta, y, X, offset = NULL) {
...@@ -402,8 +404,6 @@ V <- cov2cor(vcov(m3)) ...@@ -402,8 +404,6 @@ V <- cov2cor(vcov(m3))
corrplot.mixed(V, upper = "ellipse", col = "gray50") corrplot.mixed(V, upper = "ellipse", col = "gray50")
dev.off() dev.off()
library(plyr)
# Tamanho das covariâncias com \alpha. # Tamanho das covariâncias com \alpha.
each(sum, mean, max)(abs(V[1, -1])) each(sum, mean, max)(abs(V[1, -1]))
...@@ -421,7 +421,6 @@ crossprod(L %*% coef(m3), ...@@ -421,7 +421,6 @@ crossprod(L %*% coef(m3),
solve(L %*% vcov(m3) %*% t(L), solve(L %*% vcov(m3) %*% t(L),
L %*% coef(m3))) L %*% coef(m3)))
library(car)
# Teste de Wald para interação (poderia ser LRT, claro). # Teste de Wald para interação (poderia ser LRT, claro).
# É necessário um objeto glm. # É necessário um objeto glm.
linearHypothesis(model = m0, linearHypothesis(model = m0,
...@@ -432,9 +431,6 @@ linearHypothesis(model = m0, ...@@ -432,9 +431,6 @@ linearHypothesis(model = m0,
#----------------------------------------------------------------------- #-----------------------------------------------------------------------
# Predição com bandas de confiança. # Predição com bandas de confiança.
library(doBy)
library(multcomp)
X <- LSmatrix(m0, effect = c("umid", "K")) X <- LSmatrix(m0, effect = c("umid", "K"))
pred <- attr(X, "grid") pred <- attr(X, "grid")
...@@ -469,11 +465,9 @@ pred$pgen <- cbind(pred$pgen, ...@@ -469,11 +465,9 @@ pred$pgen <- cbind(pred$pgen,
FUN = function(x) { FUN = function(x) {
exp(pred$pgen$eta + x) exp(pred$pgen$eta + x)
})) }))
str(pred$pgen)
pred <- ldply(pred, .id = "modelo") pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo) pred <- arrange(pred, umid, K, modelo)
str(pred)
key <- list(type = "o", divide = 1, key <- list(type = "o", divide = 1,
lines = list(pch = 1:nlevels(pred$modelo), lines = list(pch = 1:nlevels(pred$modelo),
...@@ -597,14 +591,12 @@ V <- V[-1, -1] ...@@ -597,14 +591,12 @@ V <- V[-1, -1]
U <- chol(V) U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1, aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) })) FUN = function(x) { sum(x^2) }))
pred$pgen$eta <- c(X %*% coef(m3)[-1]) pred$pgen$eta <- c(X %*% coef(m3)[-1])
pred$pgen <- cbind(pred$pgen, pred$pgen <- cbind(pred$pgen,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2, apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
FUN = function(x) { FUN = function(x) {
exp(pred$pgen$eta + x) exp(pred$pgen$eta + x)
})) }))
str(pred$pgen)
pred <- ldply(pred, .id = "modelo") pred <- ldply(pred, .id = "modelo")
pred <- arrange(pred, umid, K, modelo) pred <- arrange(pred, umid, K, modelo)
...@@ -818,11 +810,9 @@ str(pred$pois) ...@@ -818,11 +810,9 @@ str(pred$pois)
# Matrix de covariância completa e sem o alpha. # Matrix de covariância completa e sem o alpha.
V <- vcov(m3) V <- vcov(m3)
V <- V[-1,-1] V <- V[-1,-1]
U <- chol(V) U <- chol(V)
aux <- sqrt(apply(X %*% t(U), MARGIN = 1, aux <- sqrt(apply(X %*% t(U), MARGIN = 1,
FUN = function(x) { sum(x^2) })) FUN = function(x) { sum(x^2) }))
pred$pgen$eta <- c(X %*% coef(m3)[-1]) pred$pgen$eta <- c(X %*% coef(m3)[-1])
pred$pgen <- cbind(pred$pgen, pred$pgen <- cbind(pred$pgen,
apply(outer(aux, qn, FUN = "*"), MARGIN = 2, apply(outer(aux, qn, FUN = "*"), MARGIN = 2,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment