-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathtestGaussian.R
53 lines (48 loc) · 1.68 KB
/
testGaussian.R
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
rm(list = ls())
setwd("~/GitHub/GraphicalModels-BayesStat")
source("Gaussian.R")
# Generate 20 decomposable graph that will be used as the true graph to generate
# 20 different datasets
number_of_trial = 5
number_of_node = 6
trueGraphs = list()
encodedList = c()
for(i in 1:number_of_trial){
while(TRUE){
graph = erdos.renyi.game(number_of_node,0.3,type="gnp",directed = FALSE)
newGraph = as_adjacency_matrix(graph, sparse = 0)
encoded = encodeGraph(newGraph)
if(isDecomposable(newGraph) & !encoded %in% encodedList){
trueGraphs[[i]] = newGraph
encodedList = c(encodedList,encoded)
break
}
}
}
mpgs = list()
maps = list()
mpg_distances = c()
map_distances = c()
count = 1
for(trueGraph in trueGraphs){
data = generateGaussianDataFromGraph(adjacencyMatrix = trueGraph, n.obs = 10000, n.variables = number_of_node)
initialCandidate = matrix(0,number_of_node,number_of_node)
print(paste("graph ", count))
chain = MetropolisHastingsGaussian(data[[2]], initialCandidate, 200, 50, 1, prior = "Binomial", p=0.3)
# Median Probability Graph
mpg = medianProbabilityGraph(chain)
mpgs[[count]] = mpg
mpg_distances = c(mpg_distances,computeSHD(trueGraph,mpg))
# Maximum a Posteriori Graph
map = maximumPosterioriGraph(chain)
maps[[count]] = map
map_distances = c(map_distances,computeSHD(trueGraph,map))
# Increase count
count = count + 1
}
x11()
par(mfrow = c(2,1))
barplot(table(mpg_distances),main = "SHD (Median Probability Graphs)", col = rainbow(length(unique(mpg_distances))))
barplot(table(map_distances),main = "SHD (Maximum a Posteriori)", col = rainbow(length(unique(map_distances))))
print(mean(mpg_distances))
print(mean(map_distances))