-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path.Rhistory
69 lines (69 loc) · 2.93 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
rm(list = ls())
rutawork = '/home/pedrohserrano/multivariate-bayesian'
install.packages('mvtnorm')
library(mvtnorm)
dde <- read.csv(paste(rutawork,'EST46114_SSVS_Data.csv',sep = ""), header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE)
rutawork = '/home/pedrohserrano/multivariate-bayesian/'
dde <- read.csv(paste(rutawork,'EST46114_SSVS_Data.csv',sep = ""), header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE)
dde$x<- (dde$x - mean(dde$x))/sqrt(var(dde$x))
# --------------------------------------------
# Analisis frecuentista
# --------------------------------------------
X <- cbind(1,dde$x,dde$z1,dde$z2,dde$z3,dde$z4,dde$z5)
Y <- dde$y
dde_mle<- glm(Y ~ -1+X, family=binomial("probit"))
dde_mle
beta_mle<- dde_mle$coef
beta0 <- rep(0,7)
Pbeta0 <- 0.25*diag(7)
beta <- rep(0,7) # Valor inicial de la cadena de Markov
n <- nrow(dde) # Numero de observaciones/individuos
z <- rep(0,n) # Valores iniciales de las variables latentes (modelo probit)
G <- 1000 # Numero de iteraciones del MCMC
# Gibbs sampler
gt <- 1
for(gt in 1:G){
eta<- X%*%beta # Predictor lineal (variable auxiliar)
z[Y==0] <- qnorm(runif(sum(1-Y), 0,pnorm(0,eta[Y==0],1)), eta[Y==0],1)
z[Y==1] <- qnorm(runif(sum(Y), pnorm(0,eta[Y==1],1),1), eta[Y==1],1)
# Muestreo de la distribucion final completa para beta
Vbeta <- solve(Pbeta0 + t(X)%*%X)
Ebeta <- Vbeta%*%(Pbeta0%*%beta0 + t(X)%*%z)
beta <- c(rmvnorm(1,Ebeta,Vbeta))
# Output
write(t(beta),file=paste(rutawork,'EST46114_SSVS_beta.out',sep = ""), ncol=7, append=T)
print(c(gt,round(beta*100)/100))
pdf(paste(rutawork,'EST46114_SSVS_beta_traces.pdf',sep = ""),width=7,height=5)
par(mfrow=c(2,1))
for(gt in 1:G){
eta<- X%*%beta # Predictor lineal (variable auxiliar)
# Muestreo de la distribucion truncada para Z
# de las distribuciones condicionales completas
z[Y==0] <- qnorm(runif(sum(1-Y), 0,pnorm(0,eta[Y==0],1)), eta[Y==0],1)
z[Y==1] <- qnorm(runif(sum(Y), pnorm(0,eta[Y==1],1),1), eta[Y==1],1)
# Muestreo de la distribucion final completa para beta
Vbeta <- solve(Pbeta0 + t(X)%*%X)
Ebeta <- Vbeta%*%(Pbeta0%*%beta0 + t(X)%*%z)
beta <- c(rmvnorm(1,Ebeta,Vbeta))
# Output
write(t(beta),file=paste(rutawork,'EST46114_SSVS_beta.out',sep = ""), ncol=7, append=T)
print(c(gt,round(beta*100)/100))
}
pdf(paste(rutawork,'EST46114_SSVS_beta_traces.pdf',sep = ""),width=7,height=5)
par(mfrow=c(2,1))
beta_out<- matrix(scan(paste(rutawork,'EST46114_SSVS_beta.out',sep = "")), ncol=7, byrow=T)
plot(beta_out[,1],type="l",xlab="iteration",
ylab="intercept (beta_1)")
plot(beta_out[,2],type="l",xlab="iteration",ylab="slope (beta_2)")
dev.off()
pdf(paste(rutawork,'EST46114_SSVS_beta_marginal_slope.pdf',sep = ""),width=7,height=5)
slp <- beta_out[101:1000,2]
par(mfrow=c(1,1))
par(mar=c(5,5,5,5))
plot(density(slp),type="l",xlab="DDE slope (beta_1)",ylab="Posterior Density",
cex=1.2)
abline(v=mean(slp))
abline(v=0.17536139,lwd=2.5) # MLE
abline(v=0.17536139 + 1.96*0.02909*c(-1,1),lwd=2.5,lty=2)
abline(v=quantile(slp,probs=c(0.025,0.975)),lty=2)
dev.off()