mc3.R 5.5 KB
Newer Older
Gilles Kratzer's avatar
Gilles Kratzer committed
1
mc3 <- function(n.var, dag.tmp, retain, ban, max.parents, sc,sc.scaled, score.cache, score, prior.choice, prior.lambda, prior.dag, verbose, scaling.factor) {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
2

Gilles Kratzer's avatar
Gilles Kratzer committed
3
4
5
6
    ## construction of neighbours list
    neighbours.list <- NULL
    alpha <- 0
    rejection <- 1
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
7

Gilles Kratzer's avatar
Gilles Kratzer committed
8
9
10
11
12
13
    ## choose mc move
    if (sum(dag.tmp) > 0) {
        if (sum(dag.tmp) < n.var * max.parents) {
            mc.move <- sample(x = c("removal", "addition", "reversal"), size = 1)
        } else {
            mc.move <- "removal"  #rep
Gilles Kratzer's avatar
Gilles Kratzer committed
14
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
15
16
    } else {
        mc.move <- "addition"  #rep
Gilles Kratzer's avatar
Gilles Kratzer committed
17
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
18

Gilles Kratzer's avatar
Gilles Kratzer committed
19
20
21
22
    switch(mc.move, removal = {
        ## Removal of an arc
        for (a in 1:n.var) {
            for (b in 1:n.var) {
Gilles Kratzer's avatar
update    
Gilles Kratzer committed
23
                if (dag.tmp[a, b] == 1 && retain[a, b]!=1) {
Gilles Kratzer's avatar
Gilles Kratzer committed
24
25
26
27
28
                  tmp <- dag.tmp
                  tmp[a, b] <- 0
                  neighbours.list[[length(neighbours.list) + 1]] <- tmp
                }
            }
Gilles Kratzer's avatar
Gilles Kratzer committed
29
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
30
    }, addition = {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
31

Gilles Kratzer's avatar
Gilles Kratzer committed
32
33
34
        ## Addition of an arc
        for (a in 1:n.var) {
            for (b in 1:n.var) {
Gilles Kratzer's avatar
update    
Gilles Kratzer committed
35
                if(((dag.tmp[a, b] == 0 & dag.tmp[b, a] == 0) & sum(dag.tmp[a, ]) <= (max.parents - 1)) && ban[a,b]!=1) {
Gilles Kratzer's avatar
Gilles Kratzer committed
36
37
38
                  tmp <- dag.tmp
                  # print(tmp)
                  tmp[a, b] <- 1
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
39
40


Gilles Kratzer's avatar
Gilles Kratzer committed
41
42
43
44
45
                  if (!identical(topoSortMAT(amat = tmp, index = FALSE), character(0))) {
                    neighbours.list[[length(neighbours.list) + 1]] <- tmp
                  }
                }
            }
Gilles Kratzer's avatar
Gilles Kratzer committed
46
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
47
    }, reversal = {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
48

Gilles Kratzer's avatar
Gilles Kratzer committed
49
50
51
        ## Reversal of an arc
        for (a in 1:n.var) {
            for (b in 1:n.var) {
Gilles Kratzer's avatar
update    
Gilles Kratzer committed
52
                if ((dag.tmp[a, b] == 1 & sum(dag.tmp[b, ]) <= (max.parents - 1)) && retain[a,b] && ban[b,a]!=1) {
Gilles Kratzer's avatar
Gilles Kratzer committed
53
54
55
56
57
58
59
60
61
62
                  tmp <- dag.tmp
                  tmp[a, b] <- 0
                  tmp[b, a] <- 1
                  if (!identical(topoSortMAT(amat = tmp, index = FALSE), character(0))) {
                    neighbours.list[[length(neighbours.list) + 1]] <- tmp
                  }
                }
            }
        }
    })
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
63
64


Gilles Kratzer's avatar
Gilles Kratzer committed
65
66
67
68
69
    ## choose of an mcmc move
    n.G <- length(neighbours.list)
    if (0 != n.G) {
        mcmc.move <- sample(x = 1:n.G, size = 1, replace = FALSE)
        dag.gprime <- neighbours.list[[mcmc.move]]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
70

Gilles Kratzer's avatar
Gilles Kratzer committed
71
72
        ## construction of neighbours list of GPRIME
        neighbours.list.gprime <- NULL
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
73

Gilles Kratzer's avatar
Gilles Kratzer committed
74
75
76
77
78
79
80
81
82
83
84
85
        switch(mc.move, removal = {
            ## Removal of an arc
            for (a in 1:n.var) {
                for (b in 1:n.var) {
                  if (dag.gprime[a, b] == 1) {
                    tmp <- dag.gprime
                    tmp[a, b] <- 0
                    neighbours.list.gprime[[length(neighbours.list.gprime) + 1]] <- tmp
                  }
                }
            }
        }, addition = {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
86

Gilles Kratzer's avatar
Gilles Kratzer committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
            ## Addition of an arc
            for (a in 1:n.var) {
                for (b in 1:n.var) {
                  if ((dag.tmp[a, b] == 0 & dag.tmp[b, a] == 0) & sum(dag.tmp[a, ]) <= (max.parents - 1)) {
                    tmp <- dag.gprime
                    # print(tmp)
                    tmp[a, b] <- 1
                    if (!identical(topoSortMAT(amat = tmp, index = FALSE), character(0))) {
                      neighbours.list.gprime[[length(neighbours.list.gprime) + 1]] <- tmp
                    }
                  }
                }
            }
        }, reversal = {
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
101

Gilles Kratzer's avatar
Gilles Kratzer committed
102
103
104
105
106
107
108
109
110
111
112
113
114
115
            ## Reveral of an arc
            for (a in 1:n.var) {
                for (b in 1:n.var) {
                  if (dag.tmp[a, b] == 1 & sum(dag.tmp[b, ]) <= (max.parents - 1)) {
                    tmp <- dag.gprime
                    tmp[a, b] <- 0
                    tmp[b, a] <- 1
                    if (!identical(topoSortMAT(amat = tmp, index = FALSE), character(0))) {
                      neighbours.list.gprime[[length(neighbours.list.gprime) + 1]] <- tmp
                    }
                  }
                }
            }
        })
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
116

Gilles Kratzer's avatar
Gilles Kratzer committed
117
        n.Gprime <- length(neighbours.list.gprime)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
118

Gilles Kratzer's avatar
Gilles Kratzer committed
119
120
121
        ## prior computation
        prior.G <- prior.Gprime <- 1
        score.G <- score.Gprime <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
122

Gilles Kratzer's avatar
Gilles Kratzer committed
123
124
125
126
127
        # user prior
        if (prior.choice == 3) {
            prior.G <- exp(-prior.lambda * sum(abs(dag.tmp - prior.dag)))
            prior.Gprime <- exp(-prior.lambda * sum(abs(dag.gprime - prior.dag)))
        }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
128

Gilles Kratzer's avatar
Gilles Kratzer committed
129
130
131
132
133
134
        for (a in 1:n.var) {
            # Koivisto priors
            if (prior.choice == 2) {
                prior.G <- prod(1/choose(n = (n.var - 1), k = sum(dag.tmp[a, ])), prior.G)
                prior.Gprime <- prod(1/choose(n = (n.var - 1), k = sum(dag.gprime[a, ])), prior.Gprime)
            }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
135

Gilles Kratzer's avatar
update    
Gilles Kratzer committed
136
            sc.tmp <- sc[score.cache$children == a,,drop = FALSE]
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
137
            score.G <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
Gilles Kratzer's avatar
Gilles Kratzer committed
138
                ])))), n.var + 1]), score.G)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
139
            score.Gprime <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[a,
Gilles Kratzer's avatar
Gilles Kratzer committed
140
141
                ])))), n.var + 1]), score.Gprime)
        }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
142
143


Gilles Kratzer's avatar
Gilles Kratzer committed
144
145
146
147
        score.Gprime.scaled <- score.dag(dag.gprime,score.cache,sc.scaled)
        score.G.scaled <- score.dag(dag.tmp,score.cache,sc.scaled)

        alpha <- min(exp(( score.Gprime.scaled - score.G.scaled) * (n.G/n.Gprime) * (prior.Gprime/prior.G)), 1)
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
148

Gilles Kratzer's avatar
Gilles Kratzer committed
149
        if(!is.numeric(alpha) | is.nan(alpha)) alpha <- 0
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
150

Gilles Kratzer's avatar
Gilles Kratzer committed
151
        score <- score.G
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
152

153
        if (!is.null(dag.gprime) && runif(1)<alpha) {
Gilles Kratzer's avatar
Gilles Kratzer committed
154
155
156
            dag.tmp <- dag.gprime
            score <- score.Gprime
            rejection <- 0
Gilles Kratzer's avatar
Gilles Kratzer committed
157
        }
Gilles Kratzer's avatar
Gilles Kratzer committed
158
    }
Gilles Kratzer's avatar
typos    
Gilles Kratzer committed
159

Gilles Kratzer's avatar
Gilles Kratzer committed
160
161
    return(list(dag.tmp = dag.tmp, score = score, alpha = alpha, rejection = rejection))
}  #EOF