Last active
August 23, 2017 06:39
-
-
Save OrenBochman/dd5ed5668cdd82e2d81c07800694469b to your computer and use it in GitHub Desktop.
Annotated Examples in R by M. Scutari and J.-B. Denis (2014)
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| --- | |
| title: '# Chapter 1: Discrete Case: Multinomial Bayesian Networks' | |
| output: | |
| html_document: default | |
| html_notebook: | |
| number_sections: yes | |
| toc: yes | |
| --- | |
| This R note book is my an annotated version of the code supplied with the first chapter of | |
| Bayesian Networks - with Examples in R by M. Scutari and J.-B. Denis (2014). It is a brief tutorial of discrete bayesian network models. The note are my own though based on my understanding of the text I have also added some extra code to show more features and better visuliza the data. However I highly recommend reading the book which covers the subject fully. | |
| ```{r} | |
| library(bnlearn) | |
| ``` | |
| # 1.1 the data | |
| About the data it is a market research survey on preference of transportation | |
| * A - age - young under 30 adult 30-60 or old over 60 | |
| * S - sex - genfer | |
| * E - education level : high school or university | |
| * O - Occupation type : self or employee | |
| * R - Residence : size of city | |
| * T - Travel prefrence : car train or other | |
| # creating structure dynamicaly | |
| # Creating the graph structure | |
| Starting with an empty graph with just a list of nodes dag stands for a directed acyclic graph | |
| ```{r} | |
| dag <- empty.graph(nodes = c("A", "S", "E", "O", "R", "T")) | |
| dag | |
| graphviz.plot(dag, shape = "circle", layout="dot") | |
| ``` | |
| Clearly a disconected set of nodes. Now add some arcs. This technique where we add one arc at a time is how graphs like social network grow dynamicaly. | |
| ```{r} | |
| dag <- set.arc(dag, from = "A", to = "E") | |
| dag <- set.arc(dag, from = "S", to = "E") | |
| dag <- set.arc(dag, from = "E", to = "O") | |
| dag <- set.arc(dag, from = "E", to = "R") | |
| dag <- set.arc(dag, from = "O", to = "T") | |
| dag <- set.arc(dag, from = "R", to = "T") | |
| dag | |
| ``` | |
| we can visulize the graph structure as follows: | |
| ```{r} | |
| graphviz.plot(dag, shape = "ellipse") | |
| ``` | |
| we can query the dag | |
| ```{r} | |
| modelstring(dag) | |
| nodes(dag) | |
| arcs(dag) | |
| ``` | |
| When we want to provide the full structure of the graph from the onset we can do so using an adjecency matrix as follows: | |
| ```{r} | |
| dag2 <- empty.graph(nodes = c("A", "S", "E", "O", "R", "T")) | |
| arc.set <- matrix(c("A", "E", | |
| "S", "E", | |
| "E", "O", | |
| "E", "R", | |
| "O", "T", | |
| "R", "T"), | |
| byrow = TRUE, ncol = 2, | |
| dimnames = list(NULL, c("from", "to"))) | |
| arcs(dag2) <- arc.set | |
| all.equal(dag, dag2) | |
| ``` | |
| # 1.3 Probabistic Representation | |
| This makes use of expert knowledge an may be called an expert system approach since we build an probablistic expert system which we can later use to query | |
| We can also aquire the structure based on factors which are a representation of the graph using a product of smaller distribution. We start by specifing the levels taken by the random variables in each factors | |
| ```{r} | |
| A.lv <- c("young", "adult", "old") | |
| S.lv <- c("M", "F") | |
| E.lv <- c("high", "uni") | |
| O.lv <- c("emp", "self") | |
| R.lv <- c("small", "big") | |
| T.lv <- c("car", "train", "other") | |
| ``` | |
| Next we supply conditional probabilities based on different factors. These encode our epert understanding of the dependencies and are typicaly elicited with the help of experts. | |
| ```{r} | |
| A.prob <- array(c(0.30, 0.50, 0.20), dim = 3, | |
| dimnames = list(A = A.lv)) | |
| A.prob | |
| S.prob <- array(c(0.60, 0.40), dim = 2, | |
| dimnames = list(S = S.lv)) | |
| S.prob | |
| O.prob <- array(c(0.96, 0.04, 0.92, 0.08), dim = c(2, 2), | |
| dimnames = list(O = O.lv, E = E.lv)) | |
| O.prob | |
| R.prob <- array(c(0.25, 0.75, 0.20, 0.80), dim = c(2, 2), | |
| dimnames = list(R = R.lv, E = E.lv)) | |
| R.prob | |
| R.prob <- matrix(c(0.25, 0.75, 0.20, 0.80), ncol = 2, | |
| dimnames = list(R = R.lv, E = E.lv)) | |
| R.prob | |
| E.prob <- array(c(0.75, 0.25, 0.72, 0.28, 0.88, 0.12, 0.64, | |
| 0.36, 0.70, 0.30, 0.90, 0.10), dim = c(2, 3, 2), | |
| dimnames = list(E = E.lv, A = A.lv, S = S.lv)) | |
| T.prob <- array(c(0.48, 0.42, 0.10, 0.56, 0.36, 0.08, 0.58, | |
| 0.24, 0.18, 0.70, 0.21, 0.09), dim = c(3, 2, 2), | |
| dimnames = list(T = T.lv, O = O.lv, R = R.lv)) | |
| ``` | |
| Now we supply our knowledge of the structure, note that this structure reflects the dependencies in the factors we used. For each andom variable we supply a cpd whose structre encodes dependence on previous variables in the dag. | |
| We later check and see that the structure is equivilent to the dynamicly created graph we introduced earlier the difference being we introduce all the structure up front. | |
| ```{r} | |
| dag3 <- model2network("[A][S][E|A:S][O|E][R|E][T|O:R]") | |
| all.equal(dag, dag3) | |
| ``` | |
| we combine the local cpd in a list called cpt and then combine it with the dag to produce a model which is essentialy an probablistic expert system. | |
| ```{r} | |
| cpt <- list(A = A.prob, S = S.prob, E = E.prob, O = O.prob, R = R.prob, T = T.prob) | |
| bn <- custom.fit(dag3, cpt) | |
| ``` | |
| we can now examine the model | |
| ```{r} | |
| nparams(bn) | |
| arcs(bn) | |
| nodes(bn) | |
| children(bn,"O") | |
| parents(bn,"O") | |
| bn$R | |
| R.cpt <- coef(bn$R) | |
| bn | |
| ``` | |
| we see that using expert knowledge allows us to model the data with just 21 parameters (sum of size of each factor) to model this network, which is far less than the 143 (product of levels of all features) required to model the overall joint distribution. | |
| we can also query additional information such as data on arcs, nodes and graph structure of diffrent nodes such as children and patents of a nodes | |
| # Learning the Joint Distribution from the data | |
| We load the data which is stored as comma seperated values. You might have to set the working directory to the location of the survey data, available from the books website. | |
| Then inspect some of it and then proceed to fit it. | |
| ```{r} | |
| survey <- read.table("survey.txt", header = TRUE) | |
| head(survey) | |
| options(digits = 3) | |
| bn.mle <- bn.fit(dag, data = survey, method = "mle") | |
| prop.table(table(survey[, c("O", "E")]), margin = 2) | |
| ################################################### | |
| ### code chunk number 32: chapter1.rnw:443-444 | |
| ################################################### | |
| bn.mle$O | |
| ################################################### | |
| ### code chunk number 33: chapter1.rnw:454-456 | |
| ################################################### | |
| bn.bayes <- bn.fit(dag, data = survey, method = "bayes", | |
| iss = 10) | |
| ################################################### | |
| ### code chunk number 34: chapter1.rnw:503-504 | |
| ################################################### | |
| bn.bayes$O | |
| ################################################### | |
| ### code chunk number 35: chapter1.rnw:526-529 | |
| ################################################### | |
| bn.bayes <- bn.fit(dag, data = survey, method = "bayes", | |
| iss = 20) | |
| bn.bayes$O | |
| ################################################### | |
| ### code chunk number 36: chapter1.rnw:629-631 | |
| ################################################### | |
| (nlevels(survey[, "T"]) - 1) * (nlevels(survey[, "E"]) - 1) * | |
| (nlevels(survey[, "O"]) * nlevels(survey[, "R"])) | |
| ################################################### | |
| ### code chunk number 37: chapter1.rnw:643-644 | |
| ################################################### | |
| ci.test("T", "E", c("O", "R"), test = "mi", data = survey) | |
| ################################################### | |
| ### code chunk number 38: chapter1.rnw:647-648 | |
| ################################################### | |
| ci.test("T", "E", c("O", "R"), test = "x2", data = survey) | |
| ################################################### | |
| ### code chunk number 39: chapter1.rnw:664-665 | |
| ################################################### | |
| ci.test("T", "O", "R", test = "x2", data = survey) | |
| ################################################### | |
| ### code chunk number 40: chapter1.rnw:669-670 | |
| ################################################### | |
| options(digits = 2) | |
| ################################################### | |
| ### code chunk number 41: chapter1.rnw:676-677 | |
| ################################################### | |
| arc.strength(dag, data = survey, criterion = "x2") | |
| ################################################### | |
| ### code chunk number 42: chapter1.rnw:734-736 | |
| ################################################### | |
| set.seed(456) | |
| options(digits = 6) | |
| ################################################### | |
| ### code chunk number 43: chapter1.rnw:738-740 | |
| ################################################### | |
| score(dag, data = survey, type = "bic") | |
| score(dag, data = survey, type = "bde", iss = 10) | |
| ################################################### | |
| ### code chunk number 44: chapter1.rnw:748-749 | |
| ################################################### | |
| score(dag, data = survey, type = "bde", iss = 1) | |
| ################################################### | |
| ### code chunk number 45: chapter1.rnw:761-764 | |
| ################################################### | |
| dag4 <- set.arc(dag, from = "E", to = "T") | |
| nparams(dag4, survey) | |
| score(dag4, data = survey, type = "bic") | |
| ################################################### | |
| ### code chunk number 46: chapter1.rnw:775-778 | |
| ################################################### | |
| rnd <- random.graph(nodes = c("A", "S", "E", "O", "R", "T")) | |
| modelstring(rnd) | |
| score(rnd, data = survey, type = "bic") | |
| ################################################### | |
| ### code chunk number 47: chapter1.rnw:790-793 | |
| ################################################### | |
| learned <- hc(survey) | |
| modelstring(learned) | |
| score(learned, data = survey, type = "bic") | |
| ################################################### | |
| ### code chunk number 48: chapter1.rnw:797-798 | |
| ################################################### | |
| learned2 <- hc(survey, score = "bde") | |
| ################################################### | |
| ### code chunk number 49: chapter1.rnw:807-808 | |
| ################################################### | |
| options(digits=3) | |
| ################################################### | |
| ### code chunk number 50: chapter1.rnw:811-812 | |
| ################################################### | |
| arc.strength(learned, data = survey, criterion = "bic") | |
| ################################################### | |
| ### code chunk number 51: chapter1.rnw:816-817 | |
| ################################################### | |
| arc.strength(dag, data = survey, criterion = "bic") | |
| ################################################### | |
| ### code chunk number 52: chapter1.rnw:824-825 | |
| ################################################### | |
| options(digits = 3) | |
| ################################################### | |
| ### code chunk number 53: chapter1.rnw:880-882 | |
| ################################################### | |
| dsep(dag, x = "S", y = "R") | |
| dsep(dag, x = "O", y = "R") | |
| ################################################### | |
| ### code chunk number 54: chapter1.rnw:890-891 | |
| ################################################### | |
| path(dag, from = "S", to = "R") | |
| ################################################### | |
| ### code chunk number 55: chapter1.rnw:895-896 | |
| ################################################### | |
| dsep(dag, x = "S", y = "R", z = "E") | |
| ################################################### | |
| ### code chunk number 56: chapter1.rnw:908-909 | |
| ################################################### | |
| dsep(dag, x = "O", y = "R", z = "E") | |
| ################################################### | |
| ### code chunk number 57: chapter1.rnw:919-921 | |
| ################################################### | |
| dsep(dag, x = "A", y = "S") | |
| dsep(dag, x = "A", y = "S", z = "E") | |
| ################################################### | |
| ### code chunk number 58: chapter1.rnw:1012-1013 | |
| ################################################### | |
| library(gRain) | |
| ################################################### | |
| ### code chunk number 59: chapter1.rnw:1019-1020 | |
| ################################################### | |
| junction <- compile(as.grain(bn)) | |
| ################################################### | |
| ### code chunk number 60: chapter1.rnw:1032-1033 | |
| ################################################### | |
| options(digits = 4) | |
| ################################################### | |
| ### code chunk number 61: chapter1.rnw:1035-1038 | |
| ################################################### | |
| querygrain(junction, nodes = "T")$T | |
| jsex <- setEvidence(junction, nodes = "S", states = "F") | |
| querygrain(jsex, nodes = "T")$T | |
| ################################################### | |
| ### code chunk number 62: chapter1.rnw:1064-1066 | |
| ################################################### | |
| jres <- setEvidence(junction, nodes = "R", states = "small") | |
| querygrain(jres, nodes = "T")$T | |
| ################################################### | |
| ### code chunk number 63: chapter1.rnw:1085-1089 | |
| ################################################### | |
| jedu <- setEvidence(junction, nodes = "E", states = "high") | |
| SxT.cpt <- querygrain(jedu, nodes = c("S", "T"), | |
| type = "joint") | |
| SxT.cpt | |
| ################################################### | |
| ### code chunk number 64: chapter1.rnw:1094-1095 | |
| ################################################### | |
| querygrain(jedu, nodes = c("S", "T"), type = "marginal") | |
| ################################################### | |
| ### code chunk number 65: chapter1.rnw:1102-1103 | |
| ################################################### | |
| querygrain(jedu, nodes = c("S", "T"), type = "conditional") | |
| ################################################### | |
| ### code chunk number 66: chapter1.rnw:1129-1130 | |
| ################################################### | |
| dsep(bn, x = "S", y = "T", z = "E") | |
| ################################################### | |
| ### code chunk number 67: chapter1.rnw:1137-1138 | |
| ################################################### | |
| SxT.ct = SxT.cpt * nrow(survey) | |
| ################################################### | |
| ### code chunk number 68: chapter1.rnw:1144-1145 | |
| ################################################### | |
| chisq.test(SxT.ct) | |
| ################################################### | |
| ### code chunk number 69: chapter1.rnw:1178-1179 | |
| ################################################### | |
| set.seed(123) | |
| ################################################### | |
| ### code chunk number 70: chapter1.rnw:1181-1183 | |
| ################################################### | |
| cpquery(bn, event = (S == "M") & (T == "car"), | |
| evidence = (E == "high")) | |
| ################################################### | |
| ### code chunk number 71: chapter1.rnw:1191-1193 | |
| ################################################### | |
| cpquery(bn, event = (S == "M") & (T == "car"), | |
| evidence = (E == "high"), n = 10^6) | |
| ################################################### | |
| ### code chunk number 72: chapter1.rnw:1206-1207 | |
| ################################################### | |
| set.seed(567) | |
| ################################################### | |
| ### code chunk number 73: chapter1.rnw:1209-1211 | |
| ################################################### | |
| cpquery(bn, event = (S == "M") & (T == "car"), | |
| evidence = list(E = "high"), method = "lw") | |
| ################################################### | |
| ### code chunk number 74: chapter1.rnw:1229-1230 | |
| ################################################### | |
| set.seed(123) | |
| ################################################### | |
| ### code chunk number 75: chapter1.rnw:1232-1234 | |
| ################################################### | |
| cpquery(bn, event = (S == "M") & (T == "car"), | |
| evidence = ((A == "young") & (E == "uni")) | (A == "adult")) | |
| ################################################### | |
| ### code chunk number 76: chapter1.rnw:1244-1247 | |
| ################################################### | |
| SxT <- cpdist(bn, nodes = c("S", "T"), | |
| evidence = (E == "high")) | |
| head(SxT) | |
| ################################################### | |
| ### code chunk number 77: chapter1.rnw:1255-1256 | |
| ################################################### | |
| options(digits = 3) | |
| ################################################### | |
| ### code chunk number 78: chapter1.rnw:1258-1259 | |
| ################################################### | |
| prop.table(table(SxT)) | |
| ################################################### | |
| ### code chunk number 79: chapter1.rnw:1295-1296 | |
| ################################################### | |
| graphviz.plot(dag) | |
| ################################################### | |
| ### code chunk number 80: chapter1.rnw:1318-1320 | |
| ################################################### | |
| hlight <- list(nodes = nodes(dag), arcs = arcs(dag), | |
| col = "grey", textCol = "grey") | |
| ################################################### | |
| ### code chunk number 81: chapter1.rnw:1325-1326 | |
| ################################################### | |
| pp <- graphviz.plot(dag, highlight = hlight) | |
| ################################################### | |
| ### code chunk number 82: chapter1.rnw:1332-1335 | |
| ################################################### | |
| edgeRenderInfo(pp) <- | |
| list(col = c("S~E" = "black", "E~R" = "black"), | |
| lwd = c("S~E" = 3, "E~R" = 3)) | |
| ################################################### | |
| ### code chunk number 83: chapter1.rnw:1347-1351 | |
| ################################################### | |
| nodeRenderInfo(pp) <- | |
| list(col = c("S" = "black", "E" = "black", "R" = "black"), | |
| textCol = c("S" = "black", "E" = "black", "R" = "black"), | |
| fill = c("E" = "grey")) | |
| ################################################### | |
| ### code chunk number 84: chapter1.rnw:1355-1356 | |
| ################################################### | |
| renderGraph(pp) | |
| ################################################### | |
| ### code chunk number 85: chapter1.rnw:1392-1394 | |
| ################################################### | |
| bn.fit.barchart(bn.mle$T, main = "Travel", | |
| xlab = "Pr(T | R,O)", ylab = "") | |
| ################################################### | |
| ### code chunk number 86: chapter1.rnw:1418-1427 | |
| ################################################### | |
| Evidence <- | |
| factor(c(rep("Unconditional",3), rep("Female", 3), | |
| rep("Small City",3)), | |
| levels = c("Unconditional", "Female", "Small City")) | |
| Travel <- factor(rep(c("car", "train", "other"), 3), | |
| levels = c("other", "train", "car")) | |
| distr <- data.frame(Evidence = Evidence, Travel = Travel, | |
| Prob = c(0.5618, 0.2808, 0.15730, 0.5620, 0.2806, | |
| 0.1573, 0.4838, 0.4170, 0.0990)) | |
| ################################################### | |
| ### code chunk number 87: chapter1.rnw:1432-1433 | |
| ################################################### | |
| head(distr) | |
| ################################################### | |
| ### code chunk number 88: chapter1.rnw:1437-1448 | |
| ################################################### | |
| barchart(Travel ~ Prob | Evidence, data = distr, | |
| layout = c(3, 1), xlab = "probability", | |
| scales = list(alternating = 1, tck = c(1, 0)), | |
| strip = strip.custom(factor.levels = | |
| c(expression(Pr(T)), | |
| expression(Pr({T} * " | " * {S == F})), | |
| expression(Pr({T} * " | " * {R == small})))), | |
| panel = function(...) { | |
| panel.barchart(...) | |
| panel.grid(h = 0, v = -1) | |
| }) | |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| show how to do inference using the expert system: | |
| * how to query a fully specified row of the joint probability | |
| * how to query using evidence |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment