########################################################################################## ## Funções de correlação na Reta ######################################################### ########################################################################################## ## Modelo exponencial exponencial <- function(phi, U){ saida <- exp(-phi*U) return(saida) } ## Modelo Gaussiano gaussiano <- function(phi, U){ saida <- exp(-phi*U^2) return(saida) } ########################################################################################## ## Verificando o comportamento da função de correlação de acordo com o valor de phi ###### ########################################################################################## require(rpanel) norm.panel <- function(panel){ curve(exponencial(x, phi = panel$phi),0,20) panel } # passar os argumentos que serão fixos, abre a janelinha panel <- rp.control() # controla a média rp.slider(panel, phi, 0, 10, initval=0, showvalue=TRUE, action=norm.panel) ## Modelo Gaussiano require(rpanel) norm.panel <- function(panel){ curve(gaussiano(x, phi = panel$phi),0,20) panel } # passar os argumentos que serão fixos, abre a janelinha panel <- rp.control() # controla a média rp.slider(panel, phi, 0, 10, initval=0, showvalue=TRUE, action=norm.panel) ########################################################################################## ### Como simular deste modelos ########################################################### ########################################################################################## tempo <- seq(0,10,l=10) mat.dist <- dist(tempo, upper=TRUE,diag=TRUE) ## Montando a matriz de variancia e covariância do modelo exponencial monta.sigma <- function(sigma, phi, mat.dist){ Sigma <- sigma^2 * exponencial(phi = phi , U = as.matrix(mat.dist)) return(Sigma) } ## Gerando um conjunto de dados require(mvtnorm) Sigma <- monta.sigma(sigma = 0.5, phi = 0.3, mat.dist = mat.dist) set.seed(123) Y <- t(rmvnorm(5, mean= rep(0,10), sigma = Sigma)) Y <- c(Y) dados <- data.frame(Y = Y, tempo = rep(tempo, 5), ID = rep(1:5, each=10)) ## Plotando require(lattice) xyplot(Y ~ tempo | ID, data=dados, type="l") ######################################################################################### ## Implementação via verossimilhança #################################################### ########################################################################################## vero.exp <- function(sigma, phi, mat.dist, dados){ dados.id <- split(dados,dados$ID) sigma <- exp(sigma) phi <- exp(phi) Sigma <- monta.sigma(sigma=sigma, phi=phi,mat.dist=mat.dist) n <- length(dados.id[[1]]$Y) m <- length(dados.id) ll <- c() for(i in 1:m){ ll[i] <- dmvnorm(dados.id[[i]]$Y, mean=rep(0,n), sigma=Sigma, log=TRUE) } saida <- -sum(ll) return(saida) } ## Usando a função dentro da mle2 require(bbmle) mat.dist <- as.matrix(dist(tempo, diag=TRUE, upper=TRUE)) vero.exp(sigma=1,phi=0.25, mat.dist=mat.dist, dados=dados) ajuste <- mle2(vero.exp, start=list(sigma=1, phi = 0.5), method="BFGS", #lower=list(sigma=1e-32, phi = 1e-32), #upper=list(sigma=Inf,phi=Inf), data=list(mat.dist=mat.dist, dados=dados)) summary(ajuste) exp(coef(ajuste)) perfil <- profile(ajuste) plot(perfil) exp(confint(perfil))