Não foi possível enviar o arquivo. Será algum problema com as permissões?
Essa é uma revisão anterior do documento!
Éder David Borges da Silva
- Graduado em Engenharia Agronômica - Eng Agronômica - UFPR
- e-mail: ederdbs@gmail.com / eder@leg.ufpr.br
Área de Interesse
Disciplinas 2011/1
Minicursos
- Análise de Experimentos de longa duração II Reunião Paranaense Ciência do Solo
Códigos (Em construção)
###-----------------------------------------------------------------###
### Reversible jump MCMC
### Modelo 1 y ~ N(b0+b1*x,sigma)
### Modelo 1 y ~ N(b0+b1*x+b2*x²,sigma)
#browseURL('http://www.icmc.usp.br/~ehlers/bayes/cap4.pdf')
# pg 76
require(MASS)#mvnorm()
require(MCMCpack)#rinvgamma() dinvgamma()
require(coda)#as.mcmc
rm(list=ls())
### Ajustar Prioris
### conferir jacobiano
rj.modelo <- function(y,x,b0,b1,sigma,b01,b11,b21,sigma1,model,mu=mu,sd=sd,mu0=mu0,
mu02=mu02,V0=V0,V02=V02,v0=v0,tau0=tau0,v02=v02,tau02=tau02){
if (model == 1){
u <- rnorm(1, mu,sd)
b0_n <- b0
b1_n <- b1
sigma_n <- sigma
b01_n <- b0 * u
b11_n <- b1 * u
b21_n <- u
sigma1_n <- sigma * (u^2)
}
if (model == 2){
u <- b21
b0_n <- b01 / u
b1_n <- b11 / u
sigma_n <- sigma1 / (u^2)
b01_n <- b01
b11_n <- b11
b21_n <- b21
sigma1_n <- sigma1
}
num <- (sum(dnorm(y,b0_n+b1_n*x,sigma_n,log=TRUE))#+
# sum(dnorm(y,mu0[1],V0[1,1],log=TRUE))+
#sum(dnorm(y,mu0[1],V0[2,2],log=TRUE))+
# sum(log(dinvgamma(y,v0,tau0)))
) * u^4
den <- (sum(dnorm(y,b01_n+b11_n*x+b21_n*x^2,sigma1_n,log=TRUE))#+
# sum(dnorm(y,mu02[1],V02[1,1],log=TRUE))+
# sum(dnorm(y,mu02[2],V02[2,2],log=TRUE))+
# sum(dnorm(y,mu02[3],V02[3,3],log=TRUE))+
# sum(log(dinvgamma(y,v02,tau02)))
) * dnorm(u,0,2)
u = runif(1, 0, 1)
if (model == 1) {
aceita = min(1, num/den)
if (u < aceita) {
model = 2
b0 <- b0_n
b1 <- b1_n
sigma <- sigma_n
}
}
if (model == 2){
aceita = min(1, den/num)
if (u < aceita) {
model = 1
b01 <- b01_n
b11 <- b11_n
b21 <- b21_n
sigma1 <- sigma1_n
}
}
if (model == 1){return(list(model = model,b0=b0,b1=b1,sigma=sigma))}
if (model == 2){return(list(model = model,b01=b01,b11=b11,b21=b21,sigma1=sigma1))}
}
rjmcmc <- function(nI, x,y,burnIN,mu=mu,sd=sd) {
chain = matrix(NA, nrow = nI, ncol = 8)
nv <- c(0,0)
chain[1,1:8] = c(1)
model = 1
n <- length(y)
###----------------------------------------------------------###
###MOdel 1
X <- model.matrix(~x)
k<-ncol(X)
#beta
mu0<-rep(0,k)
V0<-100*diag(k)
#sigma2
v0<-3
tau0<-100
#Valores iniciais
chain[1,3] <-sig2draw<- 3
invV0 <- solve(V0)
XtX <- crossprod(X,X)
Xty <- crossprod(X,y)
invV0_mu0 <- invV0 %*% mu0
###----------------------------------------------------------###
# Model 2
X2 <- cbind(1,x,x^2)
k2<-ncol(X2)
#beta
mu02<-rep(0,k2)
V02<-100*diag(k2)
#sigma2
v02<-3
tau02<-100
#Valores iniciais
chain[1,7] <- sig2draw2<- 3
invV02 <- solve(V02)
XtX2 <- crossprod(X2,X2)
Xty2 <- crossprod(X2,y)
invV0_mu02 <- invV02 %*% mu02
###----------------------------------------------------------###
for (i in 2:nI) {
if (model == 1){
# Model 1
#beta
invsig2draw <- 1/sig2draw
V1<-solve(invV0+(invsig2draw) * XtX)
mu1<-V1 %*% (invV0_mu0 + (invsig2draw)* Xty)
chain[i,1:2]<-mvrnorm(n=1,mu1,V1)
# sigma
v1<-(n+2*v0)/2
yXb <- (y-X %*% chain[i,1:2])
tyXb <-t(yXb)
tau1<-(0.5)*(tyXb %*% yXb+2*tau0)
chain[i,3] <- sig2draw <- sqrt(rinvgamma(1,v1,tau1))
}
if (model == 2){
# Model 2
#beta
invsig2draw2 <- 1/sig2draw2
V12<-solve(invV02+(invsig2draw2) * XtX2)
mu12<-V12 %*% (invV0_mu02 + (invsig2draw2)* Xty2)
chain[i,4:6]<-mvrnorm(n=1,mu12,V12)
# sigma
v12<-(n+2*v02)/2
yXb2 <- (y-X2 %*% chain[i,4:6])
tyXb2 <-t(yXb2)
tau12<-(0.5)*(tyXb2 %*% yXb2+2*tau02)
chain[i,7] <- sig2draw2 <- sqrt(rinvgamma(1,v12,tau12))
}
new <- rj.modelo(y,x,chain[i,1],chain[i,2],chain[i,3],chain[i,4],chain[i,5],chain[i,6],chain[i,7],model,mu=mu,sd=sd)
model <- new$model
if (model == 1) {
chain[i, 1] = new$b0
chain[i, 2] = new$b1
chain[i, 3] = new$sigma
nv[1] = nv[1] + 1
}
if (model == 2) {
chain[i, 4] = new$b01
chain[i, 5] = new$b11
chain[i, 6] = new$b21
chain[i, 7] = new$sigma1
nv[2] = nv[2] + 1
}
}
chain[,8] <- 1
chain[is.na(chain[,1]),8] <- 2
chain <- chain[- c(1:burnIN),]
colnames(chain) <- c('b0_1','b1_1','sigma_1','b0_2','b1_2','b2_2','sigma_2','model')
return(list(as.mcmc(na.omit(chain[,1:3])),
as.mcmc(na.omit(chain[,4:7])),
as.mcmc(na.omit(chain[,8]))))
}
x <- 1:10
y <- 10+2*x^1+rnorm(x,0,5)
plot(x,y)
res <- rjmcmc(5000,x,y,1,mu=0,sd=100)
lapply(res,summary)
plot(res[[1]])
summary(lm(y~1+I(x)))
plot(res[[2]])
summary(lm(y~1+I(x)+I(x^2)))
plot(res[[3]])
##------------------------------------------------------------------###
###-----------------------------------------------------------------###
### Regressão Beta
### pacote oficial
require(betareg)
data("FoodExpenditure", package = "betareg")
fe_beta <- betareg(I(food/income) ~ income + persons , data = FoodExpenditure)
summary(fe_beta)
###-----------------------------------------------------------------###
### log vero da regressão beta com duas covariaveis,
log.vero <- function(par,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* par[4], (1-mu)*par[4],log = TRUE))
return(ll)
}
###-----------------------------------------------------------------###
opt <- optim(c(B0=-0.5,B1=-0.51,B2=0.11,phi=35),log.vero,y=FoodExpenditure$food/FoodExpenditure$income,
x1=FoodExpenditure$income,
x2=FoodExpenditure$persons,
hessian = TRUE, control=(list(fnscale=-1)))
opt
opt$par
sqrt(-diag(solve(opt$hessian)))
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)
