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/12 17:36]
eder [section 5] 		pessoais:eder [2011/08/16 01:08]
eder [section 5]
Linha 152: Linha 152:
 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 246:
         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 257:
 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>​
  

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