mcmc.R 7.92 KB
Newer Older
Gilles Kratzer's avatar
Gilles Kratzer committed
1
################################################################################ mcmcabn.R --- Author : Gilles Kratzer Document created: 22/10/2019 -
Gilles Kratzer's avatar
Gilles Kratzer committed
2
3
4
5
6

##-------------------------------------------------------------------------
## MCMC procedure
##-------------------------------------------------------------------------

Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
7
mcmcabn <- function(score.cache = NULL, score = "mlik", data.dists = NULL, max.parents = 1, mcmc.scheme = c(100, 1000,
8
    1000), seed = 42, verbose = FALSE, start.dag = NULL, prior.dag = NULL, prior.lambda = NULL, prob.rev = 0.05, prob.mbr = 0.05, scaling = 0,
Gilles Kratzer's avatar
Gilles Kratzer committed
9
    prior.choice = 2) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
10

Gilles Kratzer's avatar
Gilles Kratzer committed
11
    #################################################### Tests
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
12
13

    .tests.mcmcabn(score.cache, data.dists, max.parents, mcmc.scheme, seed, verbose, start.dag, prior.dag, prior.lambda,
14
        prob.rev, prob.mbr, prior.choice,scaling)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
15

Gilles Kratzer's avatar
Gilles Kratzer committed
16
    ## end of tests
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
17

Gilles Kratzer's avatar
Gilles Kratzer committed
18
    ## format
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
19
    if (is.character(score))
Gilles Kratzer's avatar
Gilles Kratzer committed
20
21
        score <- tolower(score)
    score <- c("bic", "aic", "mdl", "mlik")[pmatch(score, c("bic", "aic", "mdl", "mlik"))]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
22
    if (score %!in% c("bic", "aic", "mdl", "mlik"))
Gilles Kratzer's avatar
Gilles Kratzer committed
23
24
25
        stop("A method should be provided and be one of: bic, aic, mdl or mlik.")
    if (is.matrix(prior.dag)) {
        prior.choice <- 3
Gilles Kratzer's avatar
Gilles Kratzer committed
26
    }
Gilles Kratzer's avatar
Gilles Kratzer committed
27
28
29
    if (is.matrix(prior.dag) && is.null(prior.lambda)) {
        prior.lambda <- 1
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
30

Gilles Kratzer's avatar
Gilles Kratzer committed
31
    n.var <- length(data.dists)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
32

Gilles Kratzer's avatar
Gilles Kratzer committed
33
    prob.mc3 <- 1 - (prob.rev + prob.mbr)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
34

Gilles Kratzer's avatar
Gilles Kratzer committed
35
36


Gilles Kratzer's avatar
Gilles Kratzer committed
37
38
39
40
41
42
    ## output
    out.dags <- array(data = NA, dim = c(n.var, n.var, mcmc.scheme[1] + 1))
    out.scores <- rep(NA, length = mcmc.scheme[1] + 1)
    out.alpha <- rep(NA, length = mcmc.scheme[1] + 1)
    out.method <- rep(NA, length = mcmc.scheme[1] + 1)
    out.rejection <- rep(NA, length = mcmc.scheme[1] + 1)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
43

Gilles Kratzer's avatar
Gilles Kratzer committed
44
45
    ## seeding
    set.seed(seed = seed)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
46

Gilles Kratzer's avatar
update    
Gilles Kratzer committed
47
48
49
50
51
52
53
54
55
56
    ## structural restriction

    retain <- score.cache$dag.retained
    ban <- score.cache$dag.banned

    if(is.null(retain)){retain <- matrix(data = 0, nrow = n.var, ncol = n.var)}

    if(is.null(ban)){ban <- matrix(data = 0, nrow = n.var, ncol = n.var)}


Gilles Kratzer's avatar
Gilles Kratzer committed
57
58
    ## Initializing matrix
    dag.tmp <- matrix(data = 0, nrow = n.var, ncol = n.var)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
59

Gilles Kratzer's avatar
Gilles Kratzer committed
60
    ## start zero matrix if(!is.null(start.dag)){ colnames(dag.tmp) <- rownames(dag.tmp)<-name.dag<-sample(1:n.var) }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
61

Gilles Kratzer's avatar
Gilles Kratzer committed
62
63
64
65
66
67
68
69
70
71
72
73
    ## start random matrix
    if (is.null(start.dag)) {
        start.dag <- "random"
    }
    if (is.character(start.dag) && start.dag == "random") {
        vec.tmp <- c(rep(1, max.parents), rep(0, 2 * n.var - max.parents))
        for (lines in 1:(n.var - 1)) {
            dag.tmp[1 + lines, 1:lines] <- sample(vec.tmp)[1:lines]
        }
        colnames(dag.tmp) <- rownames(dag.tmp) <- name.dag <- sample(1:n.var)
        dag.tmp <- dag.tmp[order(name.dag), order(name.dag)]
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
74

Gilles Kratzer's avatar
website    
Gilles Kratzer committed
75
76
    if (is.character(start.dag) && start.dag == "hc"){

Gilles Kratzer's avatar
Gilles Kratzer committed
77
78
      start.dag <- searchHeuristic(score.cache = score.cache,
                       score =  score,
Gilles Kratzer's avatar
website    
Gilles Kratzer committed
79
80
                       num.searches = 100,
                       max.steps = 500,
Gilles Kratzer's avatar
Gilles Kratzer committed
81
82
                       seed = seed,
                       verbose = verbose,
Gilles Kratzer's avatar
website    
Gilles Kratzer committed
83
84
85
86
87
                       start.dag = NULL,
                       dag.retained = NULL,
                       dag.banned = NULL,
                       algo = "hc",
                       tabu.memory = 10,
Gilles Kratzer's avatar
Gilles Kratzer committed
88
89
                       temperature = 0.9)

Gilles Kratzer's avatar
Gilles Kratzer committed
90
      start.dag <- start.dag[["dags"]][[which.max(x = unlist(start.dag$scores))]]
Gilles Kratzer's avatar
Gilles Kratzer committed
91

Gilles Kratzer's avatar
website    
Gilles Kratzer committed
92
93
    }

Gilles Kratzer's avatar
Gilles Kratzer committed
94
95
96
97
98
99
    if (is.matrix(start.dag)) {
        start.dag[start.dag != 0] <- 1
        diag(start.dag) <- 0
        dag.tmp <- start.dag
        colnames(dag.tmp) <- rownames(dag.tmp) <- 1:n.var
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
100

Gilles Kratzer's avatar
Gilles Kratzer committed
101

Gilles Kratzer's avatar
Gilles Kratzer committed
102
103
104
105
    ## scores
    if (score %in% c("bic", "aic", "mdl")) {
        tmp <- score.cache$node.defn[, as.numeric(colnames(dag.tmp))]
        colnames(tmp) <- as.numeric(colnames(dag.tmp))
Gilles Kratzer's avatar
Gilles Kratzer committed
106

Gilles Kratzer's avatar
Gilles Kratzer committed
107
        sc <- cbind(tmp, -score.cache[[score]])
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
108

Gilles Kratzer's avatar
Gilles Kratzer committed
109
110
    }
    if (score == "mlik") {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
111

Gilles Kratzer's avatar
Gilles Kratzer committed
112
113
        tmp <- score.cache$node.defn[, as.numeric(colnames(dag.tmp))]
        colnames(tmp) <- as.numeric(colnames(dag.tmp))
Gilles Kratzer's avatar
Gilles Kratzer committed
114

Gilles Kratzer's avatar
Gilles Kratzer committed
115
116
        sc <- cbind(tmp, score.cache[[score]])
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
117

Gilles Kratzer's avatar
Gilles Kratzer committed
118
119
    sc.scaled <- sc
    tmp.fact <- abs(max(sc[,ncol(sc)]))
120
121
    if(scaling>0){
      sc.scaled[,ncol(sc)] <- sc[,ncol(sc)]/(scaling * tmp.fact)
122
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
123

Gilles Kratzer's avatar
Gilles Kratzer committed
124
125
126
127
128
129
130
131
132

    ## scoring init
    score.init <- score.dag(dag.tmp,score.cache,sc)
    # for (a in 1:n.var) {
    #     sc.tmp <- sc[score.cache$children == as.numeric(colnames(dag.tmp)[a]), ,drop = FALSE]
    #     score.init <- sum(min(sc.tmp[which(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
    #         ])))), n.var + 1]), score.init)
    #
    # }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
133
134


Gilles Kratzer's avatar
Gilles Kratzer committed
135
136
    ## out
    out.dags[, , 1] <- dag.tmp
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
137

Gilles Kratzer's avatar
Gilles Kratzer committed
138
    out.scores[1] <- score <- score.init
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
139

Gilles Kratzer's avatar
Gilles Kratzer committed
140
    out.alpha[1] <- alpha <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
141

Gilles Kratzer's avatar
Gilles Kratzer committed
142
    out.method[1] <- "MC3"
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
143

Gilles Kratzer's avatar
Gilles Kratzer committed
144
    out.rejection[1] <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
145

Gilles Kratzer's avatar
Gilles Kratzer committed
146
    ## start mcmc search
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
147
    if (verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
148
        cat("Start MCMC burn in \n")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
149

Gilles Kratzer's avatar
Gilles Kratzer committed
150
    j <- 1
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
151

Gilles Kratzer's avatar
Gilles Kratzer committed
152
    while (j <= mcmc.scheme[3]) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
153

Gilles Kratzer's avatar
Gilles Kratzer committed
154
        ## method choice:
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
155

Gilles Kratzer's avatar
Gilles Kratzer committed
156
        method.choice <- sample(x = c("MC3", "REV", "MBR"), size = 1, prob = c(prob.mc3, prob.rev, prob.mbr))
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
157

158
159
        #if(heating=="auto"){factor <- 1/(j+1)}

Gilles Kratzer's avatar
Gilles Kratzer committed
160
        switch(method.choice, MC3 = {
Gilles Kratzer's avatar
Gilles Kratzer committed
161
            out <- mc3(n.var, (dag.tmp), retain, ban, max.parents, sc,sc.scaled, score.cache, score, prior.choice, prior.lambda, prior.dag,
Gilles Kratzer's avatar
Gilles Kratzer committed
162
163
164
165
166
167
168
                verbose)
            dag.tmp <- out$dag.tmp
            score <- out$score
        }, REV = {
            if (verbose) {
                print("REV move")
            }
Gilles Kratzer's avatar
Gilles Kratzer committed
169
            out <- REV(n.var, (dag.tmp),retain, ban, max.parents, sc,sc.scaled, score.cache, score, verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
170
171
172
173
174
175
            dag.tmp <- out$dag.tmp
            score <- out$score
        }, MBR = {
            if (verbose) {
                print("MBR move")
            }
Gilles Kratzer's avatar
Gilles Kratzer committed
176
            out <- MBR(n.var, dag.tmp,retain, ban, max.parents, sc,sc.scaled, score.cache, score, verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
177
178
            dag.tmp <- out$dag.tmp
            score <- out$score
Gilles Kratzer's avatar
updates    
Gilles Kratzer committed
179
        })
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
180

Gilles Kratzer's avatar
Gilles Kratzer committed
181
        j <- j + 1
182
    }  #EOF Burn in
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
183

Gilles Kratzer's avatar
Gilles Kratzer committed
184
185
    out.scores[1] <- score
    out.dags[, , 1] <- dag.tmp
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
186

Gilles Kratzer's avatar
Gilles Kratzer committed
187
    ### EOF: Burnin phase!
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
188
189

    if (verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
190
191
        cat("Start MCMC search \n")
    for (mcmc.search in 1:mcmc.scheme[1]) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
192
        if (verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
193
            cat("processing search ...", mcmc.search, "\n")
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
194

Gilles Kratzer's avatar
Gilles Kratzer committed
195
196
197
        ## start blind search
        j <- 0
        while (j <= mcmc.scheme[2]) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
198

Gilles Kratzer's avatar
Gilles Kratzer committed
199
            method.choice <- sample(x = c("MC3", "REV", "MBR"), size = 1, prob = c(prob.mc3, prob.rev, prob.mbr))
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
200

Gilles Kratzer's avatar
Gilles Kratzer committed
201
            switch(method.choice, MC3 = {
Gilles Kratzer's avatar
Gilles Kratzer committed
202
                out <- mc3(n.var, (dag.tmp), retain, ban, max.parents, sc,sc.scaled, score.cache, score, prior.choice, prior.lambda, prior.dag,
Gilles Kratzer's avatar
Gilles Kratzer committed
203
204
205
206
207
208
209
210
211
                  verbose)
                dag.tmp <- out$dag.tmp
                score <- out$score
                alpha <- out$alpha
                rejection <- out$rejection
            }, REV = {
                if (verbose) {
                  print("REV move")
                }
Gilles Kratzer's avatar
Gilles Kratzer committed
212
                out <- REV(n.var, (dag.tmp),retain, ban, max.parents, sc,sc.scaled, score.cache, score, verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
213
214
215
216
217
218
219
220
                dag.tmp <- out$dag.tmp
                score <- out$score
                alpha <- out$alpha
                rejection <- out$rejection
            }, MBR = {
                if (verbose) {
                  print("MBR move")
                }
Gilles Kratzer's avatar
Gilles Kratzer committed
221
                out <- MBR(n.var, (dag.tmp), retain, ban, max.parents, sc,sc.scaled, score.cache, score, verbose)
Gilles Kratzer's avatar
Gilles Kratzer committed
222
223
224
225
226
                dag.tmp <- out$dag.tmp
                score <- out$score
                alpha <- out$alpha
                rejection <- out$rejection
            })
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
227

Gilles Kratzer's avatar
Gilles Kratzer committed
228
229
            j <- j + 1
        }  #EOF mcmc search
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
230

Gilles Kratzer's avatar
Gilles Kratzer committed
231
        out.dags[, , mcmc.search + 1] <- dag.tmp
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
232

Gilles Kratzer's avatar
Gilles Kratzer committed
233
        out.scores[mcmc.search + 1] <- score
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
234

Gilles Kratzer's avatar
Gilles Kratzer committed
235
        out.alpha[mcmc.search + 1] <- alpha
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
236

Gilles Kratzer's avatar
Gilles Kratzer committed
237
        out.method[mcmc.search + 1] <- method.choice
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
238

Gilles Kratzer's avatar
Gilles Kratzer committed
239
        out.rejection[mcmc.search + 1] <- rejection
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
240

Gilles Kratzer's avatar
Gilles Kratzer committed
241
    }  #EOF mcmc search
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
242
    out <- list(dags = out.dags, scores = out.scores, alpha = out.alpha, method = out.method, rejection = out.rejection,
Gilles Kratzer's avatar
Gilles Kratzer committed
243
        iterations = mcmc.scheme[1] * (mcmc.scheme[2] + 1), thinning = mcmc.scheme[2], burnin = mcmc.scheme[3], data.dist = data.dists)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
244

Gilles Kratzer's avatar
Gilles Kratzer committed
245
    class(out) <- "mcmcabn"
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
246

Gilles Kratzer's avatar
Gilles Kratzer committed
247
    return(out)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
248

Gilles Kratzer's avatar
Gilles Kratzer committed
249
}  #eof