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

Você está aqui: start » pessoais » eder

Diferenças

Diferenças

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

Link para esta página de comparações

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/16 18:11]
eder [section 5] 		pessoais:eder [2011/09/25 20:44]
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 Genômica Ampla Via ML REML INLA
 +  * {{:​pessoais:​reml_inla.r|Script}} Modelo seleção Genótipo ambiente via REML ML 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 275: Linha 277:
 plot(grid.phi,​2*(max(opt)-opt),​type='​l'​) plot(grid.phi,​2*(max(opt)-opt),​type='​l'​)
 abline(h=3.84) 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>​
  

QR Code
QR Code pessoais:eder (generated for current page)