Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
pessoais:eder [2011/06/12 17:36] eder [section 5] |
pessoais:eder [2011/09/12 21:02] eder [section 2] |
||
---|---|---|---|
Linha 12: | Linha 12: | ||
* Estatística Espacial | * Estatística Espacial | ||
* [[http://www.leg.ufpr.br/doku.php/projetos:gem2|GEM²]] Grupo de estudos em modelos mistos | * [[http://www.leg.ufpr.br/doku.php/projetos:gem2|GEM²]] Grupo de estudos em modelos mistos | ||
+ | * {{:pessoais:inlarrblup.r|GWS}} Seleção Genomica Ampla Via ML REML INLA | ||
===== Disciplinas 2011/1 ===== | ===== Disciplinas 2011/1 ===== | ||
* [[http://www.leg.ufpr.br/doku.php/disciplinas:ce210-2010-02|CE-210: Inferência estatística II]] | * [[http://www.leg.ufpr.br/doku.php/disciplinas:ce210-2010-02|CE-210: Inferência estatística II]] | ||
Linha 152: | Linha 153: | ||
plot(x,f(x),type='l',ylim=c(0,2)) | plot(x,f(x),type='l',ylim=c(0,2)) | ||
lines(x,((2*pi)/((-D2f(0))))^0.5*f(x),col="red") | lines(x,((2*pi)/((-D2f(0))))^0.5*f(x),col="red") | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### Solução analitica, númerica e por simulação do modelo | ||
+ | # X ~ B(n,p) | ||
+ | # p ~ Beta(alfa,beta) | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | require(sfsmisc) | ||
+ | require(latticeExtra) | ||
+ | require(MASS) | ||
+ | #browseURL('http://cs.illinois.edu/class/sp10/cs598jhm/Slides/Lecture02HO.pdf') | ||
+ | |||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### grid de p | ||
+ | p <- seq(0,0.99999,by=0.001) | ||
+ | ### Priori | ||
+ | alfa <- 1 | ||
+ | beta <- 1 | ||
+ | p.priori <- dbeta(p,alfa,beta) | ||
+ | ### Verossimilhança | ||
+ | n <- 1000 | ||
+ | x <- rbinom(1,n,0.3) | ||
+ | vero <- function(p,n,x){exp(sum(dbinom(x,n,p,log=TRUE)))} | ||
+ | p.vero <- apply(matrix(p),1,vero,n=n,x=x) | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### Solução analitica | ||
+ | ### Posteriori | ||
+ | p.posteA <- dbeta(p,alfa+sum(x),beta+sum(n-x)) | ||
+ | ### Plotando | ||
+ | doubleYScale(xyplot(p.priori + p.posteA ~ p, foo, type = "l",lwd=3), | ||
+ | xyplot(p.vero ~ p, foo, type = "l",lwd=2,lty=2), | ||
+ | style1 = 0, style2 = 3, add.ylab2 = TRUE, | ||
+ | text = c("Priori", "Posteriori", "Verossimilhança"), columns = 3) | ||
+ | ### confirmando se a posteriori é uma fdp | ||
+ | integrate.xy(p,p.posteA) | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### INtegração númerica para normalização | ||
+ | ### posteriori | ||
+ | p.posteN <- (p.priori*p.vero)/(integrate.xy(p,p.priori*p.vero)) | ||
+ | ### Plotando | ||
+ | doubleYScale(xyplot(p.priori + p.posteN ~ p, foo, type = "l",lwd=2), | ||
+ | xyplot(p.vero ~ p, foo, type = "l",lwd=2,lty=2), | ||
+ | style1 = 0, style2 = 3, add.ylab2 = TRUE, | ||
+ | text = c("Priori", "Posteriori", "Verossimilhança"), columns = 3) | ||
+ | ### confirmando se a posteriori é uma fdp | ||
+ | integrate.xy(p,p.posteN) | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### Amostragem da posteriori | ||
+ | ns <- 100000 | ||
+ | theta_chapeu <- sum(x)/(n*length(x)) | ||
+ | theta_i <- rbeta(ns,alfa,beta) | ||
+ | u_i <- runif(ns,0,1) | ||
+ | crite <- u_i <= ((dbeta(theta_i,alfa,beta)*apply(matrix(theta_i),1,vero,n=n,x=x))/ | ||
+ | (dbeta(theta_chapeu,alfa,beta)*vero(theta_chapeu,n=n,x=x))) | ||
+ | a.posteriori <- theta_i[crite] | ||
+ | mean(a.posteriori,na.rm=TRUE) | ||
+ | ### Taxa Aceitação | ||
+ | sum(crite)/ns | ||
+ | ###------------------------------------------------------------### | ||
+ | ###------------------------------------------------------------### | ||
+ | ### Comparando os resultados | ||
+ | hist(a.posteriori,prob=TRUE) | ||
+ | rug(a.posteriori) | ||
+ | lines(density(a.posteriori)) | ||
+ | lines(p,p.posteA,col='red',lwd=3) | ||
+ | lines(p,p.posteN,col='blue',lty=2) | ||
+ | legend('topleft',c('Amostragem','Analitico','Númerica'),lty=c(1,1,2),col=c('black','red','blue')) | ||
+ | |||
+ | ### Intervalos via verosimilhança aproximado | ||
+ | theta_chapeu+c(-1,1)*1.96*sqrt((theta_chapeu*(1-theta_chapeu))/n) | ||
+ | ### IC amostragem | ||
+ | quantile(a.posteriori,c(0.025,0.975)) | ||
+ | ### Analitico da conjugada | ||
+ | qbeta(c(0.025,0.975),alfa+sum(x),beta+sum(n-x)) | ||
+ | ###------------------------------------------------------------### | ||
+ | ##------------------------------------------------------------### | ||
###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### | ||
### Regressão Beta | ### Regressão Beta | ||
Linha 166: | Linha 247: | ||
return(ll) | return(ll) | ||
} | } | ||
+ | |||
###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### | ||
opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income, | opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income, | ||
Linha 177: | Linha 258: | ||
summary(fe_beta) | summary(fe_beta) | ||
###-----------------------------------------------------------------### | ###-----------------------------------------------------------------### | ||
+ | log.veroP <- function(par,phi,y,x1,x2){ | ||
+ | mu <- exp((par[1] + par[2] * x1 + par[3] * x2))/(1+exp((par[1] + par[2] * x1 + par[3] * x2)))##logit^-1 | ||
+ | ll <- sum(dbeta(y, mu* phi, (1-mu)*phi,log = TRUE)) | ||
+ | return(ll) | ||
+ | } | ||
+ | |||
+ | opt <- grid.phi <- seq(20,60,l=150) | ||
+ | con <- 1 | ||
+ | for (i in grid.phi){ | ||
+ | opt[con] <- optim(c(B0=-0.5,B1=-0.51,B2=0.11),log.veroP,phi=i,y=FoodExpenditure$food/FoodExpenditure$income, | ||
+ | x1=FoodExpenditure$income, | ||
+ | x2=FoodExpenditure$persons, | ||
+ | hessian = TRUE, control=(list(fnscale=-1)))$value | ||
+ | con <- con+1 | ||
+ | } | ||
+ | |||
+ | plot(grid.phi,2*(max(opt)-opt),type='l') | ||
+ | abline(h=3.84) | ||
+ | ############################################################################################## | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ### Algumas funções para analise utilizando MCMC | ||
+ | ### topicos: | ||
+ | ### 1) Modelos com JAGS via R (Regressão linear, Logistica) | ||
+ | ### 2) dclone juntamente com JAGS | ||
+ | ### 3) MCMCsamp em modelos mistos | ||
+ | ### 4) MCMCglmm | ||
+ | ### 5) MCMCpack | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | rm(list=ls()) | ||
+ | ### JAGS | ||
+ | require(runjags) | ||
+ | ###--------------------------------------------------### | ||
+ | #### Exemplo JAGS (função run-jags) | ||
+ | ## Exemplo de regressão linear simples | ||
+ | # Simulação de dados | ||
+ | X <- 1:50 | ||
+ | Y <- rnorm(length(X), 2*X + 10, 6) | ||
+ | plot(X,Y) | ||
+ | ###--------------------------------------------------### | ||
+ | ### Inferencia por minimos quadrados | ||
+ | m0 <- lm(Y~X) | ||
+ | summary(m0) | ||
+ | ###--------------------------------------------------### | ||
+ | ### Modelo y = a*x+b | ||
+ | ### Ajustando com o JAGS | ||
+ | ### Priori para a ~ dnorm(0,0.001) | ||
+ | ### Priori para b ~ dnorm(0,0.001) | ||
+ | ### Priori para sigma ~ exp(1) | ||
+ | ### Escrevendo o modelo | ||
+ | model <- "model { | ||
+ | for(i in 1 : N){ | ||
+ | Y[i] ~ dnorm(y.est[i],sigma); | ||
+ | y.est[i] <- (a * X[i]) + b; | ||
+ | } | ||
+ | a ~ dnorm(0,0.001); | ||
+ | b ~ dnorm(0,0.001); | ||
+ | sigma ~ dexp(1); | ||
+ | }" | ||
+ | |||
+ | ### Data e valor inicial | ||
+ | data <- dump.format(list(X=X, Y=Y, N=length(X))) | ||
+ | inits0 <- dump.format(list(m=1, c=10, sigma=1)) | ||
+ | inits1 <- dump.format(list(m=0, c=0, sigma=1)) | ||
+ | # Run the model | ||
+ | m1 <- run.jags(model=model, monitor=c("a","b","sigma"), | ||
+ | data=data, inits=c(inits0,inits1),n.chains=3, plots = TRUE) | ||
+ | ### informações do objeto | ||
+ | names(m1) | ||
+ | ### Verificando a cadeia | ||
+ | plot(m1$mcmc[[1]]) | ||
+ | ### Intervalos de credibilidade | ||
+ | m1$HPD | ||
+ | ### Sumario | ||
+ | m1$summary | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ### Regressão logistica (Livro Introdução Analise Bayesiana) | ||
+ | ## pi (n/y) é a proporção de embriões com asas | ||
+ | ## ti é a variável “tempo desde a deposição dos ovos” | ||
+ | ## função de ligação logit | ||
+ | |||
+ | t <- c(5,6,8,8,10,11,16,18) | ||
+ | n <- c(34,33,33,35,30,27,33,39) | ||
+ | y <- c(6,4,23,18,28,27,33,39) | ||
+ | |||
+ | dados <- cbind(y,n-y) | ||
+ | m0 <- glm(dados~t,family=binomial(link=logit)) | ||
+ | summary(m0) | ||
+ | |||
+ | datalist <- dump.format(list(metamorfose=y,total=n,tempo=t)) | ||
+ | params <- c("beta0","beta1") | ||
+ | inicial <- dump.format(list(beta0=coef(m0)[1], | ||
+ | beta1=coef(m0)[2])) | ||
+ | ###--------------------------------------------------### | ||
+ | modmcmc <- "model{ | ||
+ | for(i in 1:length(metamorfose)){ | ||
+ | metamorfose[i] ~ dbin(p[i],total[i]) | ||
+ | logit(p[i]) <- beta0+beta1*tempo[i] | ||
+ | } | ||
+ | beta0 ~ dnorm(0,0.001) | ||
+ | beta1 ~ dnorm(0,0.001) | ||
+ | }" | ||
+ | |||
+ | ### rodando o modelo | ||
+ | modfit <- run.jags(model=modmcmc,monitor=params,data=datalist, | ||
+ | inits=inicial,n.chains=1,burnin=10000,thin=3,sample=9000,check.conv=TRUE) | ||
+ | ### Verificando a cadeia | ||
+ | plot(modfit$mcmc[[1]]) | ||
+ | ### Intervalos de credibilidade | ||
+ | modfit$HPD | ||
+ | ### Sumario | ||
+ | modfit$summary | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ### Modelos com repetição | ||
+ | x <- gl(5,4) | ||
+ | y <- rnorm(20,10,5) | ||
+ | X <- model.matrix(~x) | ||
+ | datalist <- dump.format(list(y=y,X=X)) | ||
+ | params <- c("beta1","tau") | ||
+ | inicial <- dump.format(list(beta1=rep(0,nlevels(x)),tau=1)) | ||
+ | ###--------------------------------------------------### | ||
+ | modmcmc <- "model{ | ||
+ | for(i in 1:length(y)){ | ||
+ | y[i] ~ dnorm(yb[i],tau) | ||
+ | yb[i] <- inprod(X[i,], beta1) | ||
+ | } | ||
+ | tau ~ dexp(1); | ||
+ | for(j in 1:5){ | ||
+ | beta1[j] ~ dnorm(0,0.001) | ||
+ | } | ||
+ | }" | ||
+ | ### rodando o modelo | ||
+ | modfit <- run.jags(model=modmcmc,monitor=params,data=datalist, | ||
+ | inits=inicial,check.conv=TRUE) | ||
+ | ### Verificando a cadeia | ||
+ | plot(modfit$mcmc[[1]]) | ||
+ | ### Intervalos de credibilidade | ||
+ | modfit$HPD | ||
+ | ### Sumario | ||
+ | modfit$summary | ||
+ | #### Minimos quadrados | ||
+ | summary(lm(y~x)) | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ### JAGS com Dclone | ||
+ | require(dclone) | ||
+ | ## simple regression example from the JAGS manual | ||
+ | jfun <- function() { | ||
+ | for (i in 1:N) { | ||
+ | Y[i] ~ dnorm(mu[i], tau) | ||
+ | mu[i] <- alpha + beta * (x[i] - x.bar) | ||
+ | } | ||
+ | x.bar <- mean(x[]) | ||
+ | alpha ~ dnorm(0.0, 1.0E-4) | ||
+ | beta ~ dnorm(0.0, 1.0E-4) | ||
+ | sigma <- 1.0/sqrt(tau) | ||
+ | tau ~ dgamma(1.0E-3, 1.0E-3) | ||
+ | } | ||
+ | ## data generation | ||
+ | set.seed(1234) | ||
+ | N <- 100 | ||
+ | alpha <- 1 | ||
+ | beta <- -1 | ||
+ | sigma <- 0.5 | ||
+ | x <- runif(N) | ||
+ | linpred <- model.matrix(~x) %*% c(alpha, beta) | ||
+ | Y <- rnorm(N, mean = linpred, sd = sigma) | ||
+ | ## list of data for the model | ||
+ | jdata <- list(N = N, Y = Y, x = x) | ||
+ | ## what to monitor | ||
+ | jpara <- c("alpha", "beta", "sigma") | ||
+ | ## fit the model with JAGS | ||
+ | regmod <- jags.fit(jdata, jpara, jfun, n.chains = 3) | ||
+ | ## model summary | ||
+ | summary(regmod) | ||
+ | ## data cloning | ||
+ | dcdata <- dclone(jdata, 5, multiply = "N") | ||
+ | dcmod <- jags.fit(dcdata, jpara, jfun, n.chains = 3) | ||
+ | summary(dcmod) | ||
+ | #### opções de computação paralela | ||
+ | ?jags.parfit | ||
+ | |||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ### Modelo misto - Amostrando a posteriorir via MCMC | ||
+ | require(lme4) | ||
+ | ### Resposta normal | ||
+ | m0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy) | ||
+ | summary(m0) | ||
+ | sampm0 <- mcmcsamp(m0, n = 1000) | ||
+ | HPDinterval(sampm0) | ||
+ | xyplot(sampm0) | ||
+ | qqmath(sampm0) | ||
+ | densityplot(sampm0) | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | browseURL("http://cran-r.c3sl.ufpr.br/web/packages/MCMCglmm/vignettes/CourseNotes.pdf") | ||
+ | browseURL("http://cran-r.c3sl.ufpr.br/web/packages/MCMCglmm/vignettes/Overview.pdf") | ||
+ | browseURL("http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.160.5098&rep=rep1&type=pdf") | ||
+ | require(MCMCglmm) | ||
+ | ?MCMCglmm | ||
+ | data(PlodiaPO) | ||
+ | ### Modelo normal | ||
+ | model1<-MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=FALSE) | ||
+ | summary(model1) | ||
+ | plot.MCMCglmm(model1) | ||
+ | ### MOdelo binomial | ||
+ | model2 <- MCMCglmm(cbind(Pupated, Infected) ~ 1,random=~FSfamily, family = "multinomial2", | ||
+ | data = PlodiaR, verbose = FALSE) | ||
+ | plot.MCMCglmm(model2) | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | require(MCMCpack) | ||
+ | ### Modelo Poisson | ||
+ | counts <- c(18,17,15,20,10,20,25,13,12) | ||
+ | outcome <- gl(3,1,9) | ||
+ | treatment <- gl(3,3) | ||
+ | posterior <- MCMCpoisson(counts ~ outcome + treatment) | ||
+ | plot(posterior) | ||
+ | summary(posterior) | ||
+ | ###--------------------------------------------------### | ||
+ | ###--------------------------------------------------### | ||
+ | ############################################################################################## | ||
</code> | </code> | ||