bayesian-learning-with-R/missingExOpenbugs.R at master · NCSU-STATS/bayesian-learning-with-R · GitHub

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
setwd("d:/brucebcampbell-git/bayesian-learning-with-R")
load("E2.RData")

x<- as.matrix(Y1)

stacks_dat <- list(x=x,p = 6,   N = 365)

all_data           <- cbind(x)


mlr_model <- function(){

    # Missing data model for x
  for(i in 1:N){
    x[i,1:p]~dmnorm(x_mn[],x_prec[,])
  }

  # Priors for missing-data model parameters
  for(j in 1:p){
    x_mn[j]~dnorm(0,0.01)
  }
  x_prec[1:p,1:p]~dwish(R[,],k)
  x_cov[1:p,1:p]<-inverse(x_prec[,])

  k<-p+0.1
  for(j1 in 1:p){for(j2 in 1:p){R[j1,j2]<-0.1*equals(j1,j2)}} #R is diagonal
}

library(rjags)
library(coda)

library(R2OpenBUGS)

mlr_inits <- function() {
  list(intercept = mean(stacks_dat$Y,na.rm=T), beta = rnorm(stacks_dat$p), tau = 0.01)
}

samps <- bugs(data = stacks_dat,
              inits = mlr_inits,
              parameters.to.save = c("x_cov", "x"),
              model.file = mlr_model, debug=FALSE,
              n.chains = 1, n.burnin=1000, n.iter = 10000, n.thin=2, DIC=F)

library(coda)

out.coda <- read.bugs(samps)
#save(out.coda,file = "out.coda_OPENBUGS_VAR1.RData")

if(n.chains > 1)
{
  g <- matrix(NA, nrow=nvar(out.coda), ncol=2)
  for (v in 1:nvar(out.coda)) {
    g[v,] <- gelman.diag(out.coda[,v])$psrf
  }
  count.coeff.gt <- sum(g[,1]>1.1)
  count.coeff.gt
  plot(g[,1],main="Gelman-Rubin ")
}

chains.ess <- lapply(out.coda,effectiveSize)

first.chain.ess <- chains.ess[1]
plot(unlist(first.chain.ess), main="Effective Sample Size")

chain <- out.coda[[1]]

Y1.posterior.modes <- list()

for( i in 1:(365*6+36) )
{
  colname <- colnames(chain)[i]
  if(grepl("Sigma",colname)  )
  {
    print(paste("skipping ",colname,sep=""))
    next
  }
  samples <- chain[,i]
  samplesId <- samples<4
  red <- samples[samplesId]
  Y1.posterior.modes[colname] <-mlv(red)$M
}

theta.map <-  matrix(unlist(Y1.posterior.modes),ncol=6, byrow=FALSE)

unscaled.theta.map <- ( theta.map +colMeans(Y1,na.rm = TRUE)) *  apply(Y1, 2,sd,na.rm = TRUE)

image(Y1, main="Y1")

image(unscaled.theta.map, main="Imputed")

SigmaVAR.posterior.modes <- list()

for( i in 1:(365*6+36) )
{
  colname <- colnames(chain)[i]
  if(grepl("x_cov",colname) ==FALSE )
  {
    next
  }
  samples <- chain[,i]
  SigmaVAR.posterior.modes[colname] <-mlv(samples)$M
}

Sigma.map <-  matrix(unlist(SigmaVAR.posterior.modes),ncol=6, byrow=FALSE)
image(Sigma.map, main = "MAP Sigma")


b <- unscaled.theta.map[!is.na(Y1)]
c <- Y1[!is.na(Y1)]

plot(b,c)