Maio 2012 - R-br - DINF UFPR

Dúvida em erro na função Filtro de Kalman
by Victor Eduardo 14 Mai '12

14 Mai '12

Galera, estou com um problema que não estou entendendo o motivo do erro que está na função filtro_kalman. Estou realizando um exercício que é para construir uma função que realiza o filtro de kalman. Utilizei o modelo de tendencia linear local no filtro de kalman: #Função para gerar as matrizes do sistema para o #Modelo de tendencia linear local matriz_Tll<-function(fi,sig2_eta,sig2_eps,sig2_qui,n){ #Definindo as dimensões das matrizes do sistema ar_Z<-array( ,dim=c(1,2,n)) ar_T<-array( ,dim=c(2,2,n)) ar_C<-array( ,dim=c(2,1,n)) ar_D<-array( ,dim=c(1,1,n)) ar_R<-array( ,dim=c(2,2,n)) ar_H<-array( ,dim=c(1,1,n)) ar_Q<-array( ,dim=c(2,2,n)) #Iniciando um 'for' para 'construir' as matrizes do sistema for(t in 1:n){ ar_Z[ , ,t]<-matrix(c(1,0),nc=2) ar_T[ , ,t]<-matrix(c(1,0,1,fi),nc=2) ar_C[ , ,t]<-matrix(c(0,0)) ar_D[ , ,t]<-matrix(0) ar_R[ , ,t]<-matrix(c(1,0,0,1),nc=2) ar_H[ , ,t]<-matrix(sig2_eps,nc=1) ar_Q[ , ,t]<-matrix(c(sig2_eta,0,0,sig2_qui),nc=2) } return(list(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q)) } fi<-0.5 #com |fi|<1 temos um modelo de tendencia amortecida sig2_eta<-4 sig2_eps<-6 sig2_qui<-10 n<-100 #número pequeno apenas para teste mat_sist<-matriz_Tll(fi,sig2_eta,sig2_eps,sig2_qui,n) ar_Z<-mat_sist[[1]] ar_T<-mat_sist[[2]] ar_C<-mat_sist[[3]] ar_D<-mat_sist[[4]] ar_R<-mat_sist[[5]] ar_H<-mat_sist[[6]] ar_Q<-mat_sist[[7]] #Função para modelo em espaço de estado em geral gerador_EE<-function(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1){ n<-dim(ar_Z)[3] #nº de matrizes do array Z ou D ou H p<-dim(ar_H)[1] #nº linhas de H, D ou Z m<-dim(ar_T)[1] #nº linhas de T, C ou R r<-dim(ar_Q)[1] #nº linhas de Q ar_Y<-array(,dim=c(p,1,n)) eps<-array(,dim=c(p,1,n)) eta<-array(,dim=c(r,1,n)) alfa<-array(,dim=c(m,1,n+1)) alfa[,,1]<-a1+chol(p1)%*%matrix(rnorm(m,0,1),nc=1) for(t in 1:n){ eps[,,t]<-0+chol(ar_H[,,t])%*%matrix(rnorm(p,0,1),nc=1) ar_Y[,,t]<-ar_Z[,,t]%*%alfa[,,t]+ar_D[,,t]+eps[,,t] eta[,,t]<-0+chol(ar_Q[,,t])%*%matrix(rnorm(r,0,1),nc=1) alfa[,,t+1]<-ar_T[,,t]%*%alfa[,,t]+ar_C[,,t]+ar_R[,,t]%*%eta[,,t] } return(ar_Y) } a1<-matrix(1:2,nc=1) p1<-matrix(c(3,1,1,2),nc=2) geracao<-gerador_EE(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1) #Gráfico dos valores gerados para Y plot(geracao,type='l') #Filtro de Kalman install.packages('MASS') require(MASS) Filtro_Kalman<-function(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1){ n<-dim(ar_Z)[3] #nº de matrizes do array Z ou D ou H p<-dim(ar_H)[1] #nº linhas de H, D ou Z m<-dim(ar_T)[1] #nº linhas de T, C ou R r<-dim(ar_Q)[1] #nº linhas de Q ar_Y<-array(,dim=c(p,1,n)) ar_U<-array(,dim=c(p,1,n)) ar_F<-array(,dim=c(p,p,n)) ar_AA<-array(,dim=c(m,1,n)) ar_PP<-array(,dim=c(m,m,n)) ar_A<-array(,dim=c(m,1,n+1)) ar_P<-array(,dim=c(m,m,n+1)) ar_A[,,1]<-a1 ar_P[,,1]<-p1 ar_Y<-gerador_EE(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1) for(t in 1:2){ ar_U[,,t]<-ar_Y[,,t]-ar_Z[,,t,drop=FALSE]%*%ar_A[,,t,drop=FALSE]-ar_D[,,t,drop=FALSE] ar_F[,,t]<-ar_Z[,,t,drop=FALSE]%*%ar_P[,,t,drop=FALSE]%*%t(ar_Z[,,t,drop=FALSE])+ar_H[,,t,drop=FALSE] ar_AA[,,t]<-ar_A[,,t,drop=FALSE]+ar_P[,,t,drop=FALSE]%*%ar_Z[,,t,drop=FALSE]%*%ginv(ar_F[,,t,drop=FALSE])%*%ar_U[,,t,drop=FALSE] ar_PP[,,t]<-ar_P[,,t,drop=FALSE]-ar_P[,,t,drop=FALSE]%*%t(ar_Z[,,t,drop=FALSE])%*%ginv(ar_F[,,t,drop=FALSE])%*%ar_Z[,,t,drop=FALSE]%*%ar_P[,,t] ar_A[,,t+1]<-ar_T[,,t,drop=FALSE]%*%ar_AA[,,t,drop=FALSE]+ar_C[,,t,drop=FALSE] ar_P[,,t+1]<-ar_T[,,t,drop=FALSE]%*%ar_PP[,,t,drop=FALSE]%*%t(ar_PP[,,t,drop=FALSE])+ar_R[,,t,drop=FALSE]%*%ar_Q[,,t,drop=FALSE]%*%t(ar_R[,,t,drop=FALSE]) } return(list(ar_A,ar_P)) } Filtro_Kalman(ar_Z,ar_T,ar_C,ar_D,ar_R,ar_H,ar_Q,a1,p1)

1 0

Correlação entre vetores de comprimentos diferentes
by Tito Conte 13 Mai '12

13 Mai '12

Senhores, Preciso saber quais posições de um vetor de 10000 elementos são maiores do que os valores de um vetor de 203900 elementos (há NANs e elementos repetidos) e correlacionar com um vetor resultado do mesmo processo (comparar latitudes e longitudes, de duas tabelas distintas e observar em que posição a condição é igual) o problema é que não consigo fazer essa comparação pois um vetor é menor que o outro, há alguma forma de forçar essa comparação ou a única solução eh um loop gigantesco de 203900 steps? Tito Conte

3 5

Dúvida Plot
by Sérgio Henrique almeida da silva ju 12 Mai '12

12 Mai '12

Olá Como posso plotar uma série em um mesmo gráfico? Ou seja, eu gostaria de plotar a série de PE, AL e SP em um mesmo gráfico. A unica coisa que consegui até agora foi: PE.ts <- ts(dados[,2], start = 2001, end=2009, freq = 1) AL.ts <- ts(dados[,3], start = 2001, end=2009, freq = 1) SP.ts <- ts(dados[,4], start = 2001, end=2009, freq = 1) layout(1:3) plot(PE.ts,col="green",lwd=3,xlab="Ano") plot(AL.ts,col="green",lwd=3,xlab="Ano") plot(SP.ts,col="green",lwd=3,xlab="Ano") Porém eu gostaria de que as linhas estivessem em um mesmo gráfico > dput(dados) structure(list(ano = c(2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L), PE = c(201L, 191L, 240L, 271L, 275L, 289L, 284L, 288L, 254L, 8735L, 9944L, 9533L, 4570L, 5106L, 4461L, 4107L, 2386L, 1114L), AL = c(216L, 260L, 62L, 56L, 71L, 60L, 82L, 61L, 69L, 497L, 1595L, 349L, 282L, 268L, 300L, 166L, 104L, 33L), SP = c(144L, 107L, 113L, 142L, 112L, 124L, 114L, 112L, 125L, 4189L, 2787L, 2738L, 2247L, 1778L, 1144L, 1549L, 1197L, 1147L), tipo = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("SIM", "SINAN"), class = "factor")), .Names = c("ano", "PE", "AL", "SP", "tipo"), class = "data.frame", row.names = c(NA, -18L)) Abraços -- Sérgio Henrique Almeida da Silva Junior Doutorando em Epidemiologia em Saúde Pública Escola Nacional de Saúde Pública Sérgio Arouca - ENSP/FIOCRUZ http://lattes.cnpq.br/1611345552843383 Tel: (21) 94429486/78101651 id: 123*20942

5 5

Tukey para parcelas subdividad
by Gilson Sanchez 12 Mai '12

12 Mai '12

Muito obrigado pessoal pelas dicas fiquei interessado no pacote ExpDes, mas da pra ajudar no seguinte mensajem ------------------------------------------------------------------------ Legend: FACTOR 1 (plot): solo FACTOR 2 (split-plot): tempo ------------------------------------------------------------------------ Erro em `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels * * * * *O Script que estou tentando utilizar é:* split2.rbd(dados$solo, dados$tempo, dados$bloco, dados$a, quali=c(TRUE, FALSE), mcomp = "tukey", fac.names = c("solo", "tempo"), sigT = 0.05, sigF = 0.05) Aguardo resposta -- Bach. Agr.Gilson Sánchez Chia, MSc. Laboratório de Fisiologia Vegetal Embrapa Amazônia Ocidental Fone: (92) 3303-7841 / 3303-7919 http://lattes.cnpq.br/3684553005529634 --- Antes de imprimir, pense na sua responsabilidade com o Meio Ambiente. Please consider the environment before printing this email. Renew, Reduce and Recycle. The Planet Knows.*** Salva un Árbol.* * *** * PIENSA**** EN EL MUNDO PARA TUS HIJOS** *

4 6

Modelo Bayesiano
by Sérgio Henrique almeida da silva ju 11 Mai '12

11 Mai '12

Olá Estou tentando estimar a prevalência de uma doença (True Prevalence), através de um modelo Bayesiano (Joseph L, Gyorkos TW, Coupal L Bayesian estimation of disease prevalence and the parameters of diagnostic tests in the absence of a gold standard Am J Epidemiol 1995,141:263-72.) Para isso estou usando a biblioteca (R2WinBUGS) e modelando através de distribuição Binomial Negativa ##Parameterized by p and r: model { for( i in 1 : N ) { y[i] ~ dnegbin( p , r ) } p ~ dbeta(1.001,1.001) r ~ dgamma(0.01,0.01) m <- r*(1-p)/p v <- r*(1-p)/(p*p) } Porém quando rodo o modelo: N <-nrow(dados1) y <-dados1$esqui data <- list("N", "y") inits <- function(){list(p=rbeta(N,1.001,1.001),r=rgamma(1,0.01,0.01))} esqui.sim <- bugs(data, inits, model.file = "esqui.bugs", parameters = c("p", "r"), n.chains = 2, n.iter = 1000, bugs.directory = "C:/Users/sergio/Downloads/winbugs14",debug=T) Aparece o seguinte erro no log: display(log) check(C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/esqui.bugs.txt) model is syntactically correct data(C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/data.txt) expected collection operator c compile(2) inits(1,C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/inits1.txt) command #Bugs:inits cannot be executed (is greyed out) inits(2,C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/inits2.txt) command #Bugs:inits cannot be executed (is greyed out) gen.inits() command #Bugs:gen.inits cannot be executed (is greyed out) thin.updater(1) update(500) command #Bugs:update cannot be executed (is greyed out) set(p) command #Bugs:set cannot be executed (is greyed out) set(r) command #Bugs:set cannot be executed (is greyed out) set(deviance) command #Bugs:set cannot be executed (is greyed out) dic.set() command #Bugs:dic.set cannot be executed (is greyed out) update(500) command #Bugs:update cannot be executed (is greyed out) coda(*,C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/coda) command #Bugs:coda cannot be executed (is greyed out) stats(*) command #Bugs:stats cannot be executed (is greyed out) dic.stats() DIC history(*,C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/history.odc) command #Bugs:history cannot be executed (is greyed out) save(C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/log.odc) save(C:/Users/sergio/AppData/Local/Temp/Rtmp4SGBkz/log.txt) O que estou fazendo de errado? Abraços -- Sérgio Henrique Almeida da Silva Junior Doutorando em Epidemiologia em Saúde Pública Escola Nacional de Saúde Pública Sérgio Arouca - ENSP/FIOCRUZ http://lattes.cnpq.br/1611345552843383 Tel: (21) 94429486/78101651 id: 123*20942

3 3

Ler dados do epi info
by Emerson Cotta Bodevan 11 Mai '12

11 Mai '12

Prezados, como faço para importar dados do epi info para o R.? -- *Emerson*

4 4

Simulação de um DQL no R
by Edmundo Caetano 11 Mai '12

11 Mai '12

Prezados, Pretendo simular um experimento em quadrados latinos no R. Como devo proceder??? Edmundo Caetano

4 8

Problemas com comando witch
by Tito Conte 11 Mai '12

11 Mai '12

tenho dois vetores: p=c(2,3,4,5,6,7,8,9) r=c(3,4,5) quando uso which(p==r) ou which(r==p) ele me retorna as posições onde r ocorre em p como fazer para saber onde p ocorre em r? obrigado (já dei um google nisso e não encontrei ) Tito Conte

4 7

Escalas
by Eric Keven 11 Mai '12

11 Mai '12

Boa noite, Como eu faço para alterar a escala dos gráficos e fazer com que o programacomece num valor que eu fixar para a escala?Obrigado! Att, Eric Keven Silva Coordenador do Núcleo de Estudos de Processos da Indústria de Alimentos (PROALI) Graduando em Engenharia de Alimentos Universidade Federal de Lavras

5 4

Comparações múltiplas usando MMC (mean–mean multiple comparisons) plots
by ASANTOS 10 Mai '12

10 Mai '12

*Bom dia pessoal, Estou tendo problemas em visualizar os resultados de uma comparação múltipla usando GLM, com o uso MMC plots do pacote HH. Bom, tenho quatro grupos funcionais de cupins que são humívoros (humi), ceifadores(ceifa), intermediários (inter) e xilófagos (xilo) que são meus tratamentos e como variável resposta o número de ocorrência de espécies de cada grupo (acorr). Para coleta da variável resposta, fizemos 6 transectos que são as repetições. Para ver se existe diferença entre os grupos funcionais fiz:* # mod2<-glm(ocorr~grupo, family="quasipoisson", data=trat1_90antes)## Um plot(mod2) prévio mostrou que quasipoisson foi a distribuição mais adequada aos dados que poisson summary(mod2) glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Deviance Residuals: Min 1Q Median 3Q Max -1.9148 -0.4946 0.1213 0.4178 1.3813 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.5041 0.1624 9.260 1.13e-08 *** grupohumi -1.0986 0.3249 -3.382 0.00296 ** grupointer -0.8979 0.3019 -2.974 0.00750 ** grupoxilo -0.6568 0.2780 -2.363 0.02838 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 0.7123377) Null deviance: 28.663 on 23 degrees of freedom Residual deviance: 16.889 on 20 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 anova(mod2, test="Chi") ## Anova do modelo Analysis of Deviance Table Model: quasipoisson, link: log Response: ocorr Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 23 28.663 grupo 3 11.774 20 16.889 0.0008834 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 plot(mod2) ##Teste de comparação múltipla summary(glht(mod2, linfct=mcp(grupo="Tukey"))) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Linear Hypotheses: Estimate Std. Error z value Pr(>|z|) humi - ceifa == 0 -1.0986 0.3249 -3.382 0.00405 ** inter - ceifa == 0 -0.8979 0.3019 -2.974 0.01491 * xilo - ceifa == 0 -0.6568 0.2780 -2.363 0.08255 . inter - humi == 0 0.2007 0.3794 0.529 0.95122 xilo - humi == 0 0.4418 0.3606 1.225 0.60632 xilo - inter == 0 0.2412 0.3401 0.709 0.89176 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method) > ##Resultados > tapply(trat1_90antes$ocorr, trat1_90antes$grupo,mean, na.rm = TRUE)##Media dos tratamentos ceifa humi inter xilo 4.500000 1.500000 1.833333 2.333333 > tapply(trat1_90antes$ocorr,trat1_90antes$grupo,sd,na.rm=TRUE)/sqrt(tapply(trat1_90antes$ocorr,trat1_90antes$grupo,length))###Erros padrão da media ceifa humi inter xilo 0.6708204 0.5000000 0.5426274 0.3333333 *Aqui tornou-se difícil saber a diferença entre os grupos funcionais em termos de atribui a, b, c, as médias, então resolvi fazer uma análise gráfica usando o pacote HH:* > ##Visualizando as comparações > contrast.1<- + if.R(r=glht(mod2, linfct=mcp(grupo="Tukey")), + s=multicomp(mod2, plot=FALSE)) > plot(contrast.1) *Que me deu um gráfico pouco informativo, porém seguindo o help do mmc {HH}, encontrei que o gráfico de comparações deveria ser feito usando as funções glht.mmc e multicomp.mmc, que requerem o uso de aov(), então fiz:* > > ## Contraste por mmc > contrast.aov <- aov(ocorr ~ grupo, data=trat1_90antes) > contrast.mmc <- + if.R(r=glht.mmc(contrast.aov, linfct=mcp(grupo="Tukey")), + s=multicomp.mmc(contrast.aov, plot=FALSE)) > plot(contrast.mmc) > # *Pois bem, agora tenho a comparação entre os grupos, mas não utilizando a distribuição de erros quasipoisson e não sei o que fazer, pois o pacote HH não me permite utilizar os objetos criados usando glm(), para se fazer as comparações múltiplas. Alguém do grupo teria alguma sugestão, Obrigado, Alexandre* -- Alexandre dos Santos Engenheiro Florestal, Dr. Universidade Federal de Lavras Departamento de Entomologia Laboratório de Entomologia Florestal Caixa Postal 3037 37200-000 - Lavras/MG Fone: +55 (35) 9223-0304

2 2