From 5a3b633ab446ef153ca4c62c9a71d3e3cdb52e73 Mon Sep 17 00:00:00 2001 From: Kostis Gourgoulias Date: Sat, 8 Jun 2019 12:54:11 +0100 Subject: [PATCH] It runs now but results are still weird --- bayes_net.py | 106 ++++++++++++++++++++++++++++++- update_proposals.py | 150 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 253 insertions(+), 3 deletions(-) create mode 100644 update_proposals.py diff --git a/bayes_net.py b/bayes_net.py index 1d532dc..95853d4 100644 --- a/bayes_net.py +++ b/bayes_net.py @@ -2,11 +2,111 @@ Classes to represent Bayesian networks. """ from toposort import toposort_flatten as flatten +from misc import dict_to_string +from scipy import rand, prod + + +class BayesNet(object): + """ + Object to hold and evaluate probabilities + for BayesNets described by cpts. + """ + + def __init__(self, graph, cpt): + """ + graph: dictionary of form {child:{parent1, parent2},}. + Expects {None} for the parents of root nodes. + cpt: dictionary, holds conditional probabilities for the + nodes of the network. In general, the form is expected to be: + cpt = {child: prior_probs, child:{parent+value: probs, + parent+value: probs}} + Example: + cpt = {"A":[0.2,0.8], B:{"A0":..., "A1":...}} + The values next to the node name correspond to the values + of the parent node. + """ + + self.nodes = flatten(graph) + self.nodes = self.nodes[1::] + self.graph = graph + self.cpt = cpt + + def joint_prob(self, node_values): + """ + Calculate the joint probability of an instantiation + of the graph. + Input: + node_values: dictionary, assumed to be of form {node1:value1, + node2:value2, ...} + """ + + result = 1.0 + for node in node_values: + if self.is_root_node(node): + # root node + result *= self.prior(node, node_values[node]) + else: + result *= self.cond_prob(node, node_values[node], node_values) + + return result + + def cond_prob(self, child, state, all_vals): + """ + Evaluates the conditional probability + P(child = state | all_vals) + by looking up the values from the Icpt table. + """ + parents = {key: int(all_vals[key]) for key in self.graph[child]} + key = dict_to_string(parents) + result = self.cpt[child][key][state] + return result + + def prior(self, node, value): + """ + Returns the prior of a root node. + """ + + result = None + if self.is_root_node(node): + result = self.cpt[node][value] + + return result + + def sample(self, set_nodes={}): + """ + Generate single sample from BN. + This only assumes binary variables. + """ + + # sample all but the already set nodes + nodes = (n for n in self.nodes if n not in set_nodes) + + sample = set_nodes.copy() + + for node in nodes: + if self.is_root_node(node): + p = self.prior(node, True) + else: + p = self.cond_prob(node, True, sample) + + sample[node] = int(rand() < p) + return sample + + def msample(self, num_of_samples=100): + """ + Generate multiple samples. + """ + + samples = [None] * num_of_samples + for i in range(num_of_samples): + samples[i] = self.sample() + return samples + + def is_root_node(self, node): + result = (self.graph[node] == {None}) + return result -""" -NOTE We can use this to represent a bayesian network! -""" class BNNoisyORLeaky(BayesNet): """ diff --git a/update_proposals.py b/update_proposals.py new file mode 100644 index 0000000..e92f399 --- /dev/null +++ b/update_proposals.py @@ -0,0 +1,150 @@ +from misc import weight_average, string_to_dict, char_fun + +from scipy import log, exp + + +def update_proposal_cpt(proposal, samples, weights, index, graph, + evidence_parents, eta_rate): + """ + Updates current proposal given the new data. + + Arguments + ========= + samples: the current samples to use. + index: the current index (used to pick the weight) + """ + + # Initialize weighted estimator + wei_est = weight_average(samples, weights) + + # estimate CPT table from samples + for node in evidence_parents: + if node is None: + continue + elif proposal.is_root_node(node): + # root node + def f(sample): + res1 = char_fun(sample, {node: 0}) + return res1 + + p, _ = wei_est.eval(f) + + proposal.cpt[node][0] += eta_rate(index) * ( + p - proposal.cpt[node][0]) + proposal.cpt[node][1] += eta_rate(index) * ( + 1 - p - proposal.cpt[node][1]) + else: + # rest of the nodes + for key in proposal.cpt[node]: + + parent_dict = string_to_dict(key) + + def f(sample): + res1 = char_fun(sample, {node: 0}) + res2 = char_fun(sample, parent_dict) + return res1 * res2 + + def g(sample): + res2 = char_fun(sample, parent_dict) + return res2 + + p, _ = wei_est.eval(f) + q, _ = wei_est.eval(g) + + if abs(p - q) < 1e-10: + ratio = 1 + else: + ratio = p / q + + proposal.cpt[node][key][0] += eta_rate(index) * ( + ratio - proposal.cpt[node][key][0]) + proposal.cpt[node][key][1] += eta_rate(index) * ( + 1 - ratio - proposal.cpt[node][key][1]) + + return proposal + + +def update_proposal_lambdas(proposal, samples, weights, index, graph, + evidence_parents, eta_rate): + """ + Updates current proposal given the new data. + + Arguments + ========= + samples: the current samples to use. + index: the current index (used to pick the weight) + """ + + # Initialize weighted estimator + wei_est = weight_average(samples, weights) + + # estimate CPT table from samples + for child in evidence_parents: + + if child is None: + continue + elif proposal.is_root_node(child): + # root child -- update priors using current samples + + def f(sample): + res1 = char_fun(sample, {child: 1}) + return res1 + p, _ = wei_est.eval(f) + + proposal.prior_dict[child][0] += eta_rate(index) * ( + 1 - p - proposal.prior_dict[child][0]) + proposal.prior_dict[child][1] += eta_rate(index) * ( + p - proposal.prior_dict[child][1]) + + else: + # rest of the childs -- lambdas + parents = [ident for ident in graph[child]] + for parent in proposal.lambdas[child]: + state_vec = {p: False for p in parents} + if parent == "leak_node": + def f(sample): + return char_fun(sample, state_vec) + + q, _ = wei_est.eval(f) + + state_vec[child] = False + + def f(sample): + return char_fun(sample, state_vec) + + p, _ = wei_est.eval(f) + if abs(p) < 1e-16 or abs(q) < 1e-16: + # Do not update if probabilities are + # too small + continue + ratio = exp(log(p) - log(q)) + + proposal.lambdas[child]["leak_node"] += eta_rate(index) * ( + ratio - proposal.lambdas[child]["leak_node"]) + else: + # TODO: something is wrong with this part + # and it doesn't estimate correctly + state_vec[parent] = True + + def f(sample): + return char_fun(sample, state_vec) + + q, _ = wei_est.eval(f) + state_vec[child] = False + + def f(sample): + return char_fun(sample, state_vec) + + p, _ = wei_est.eval(f) + + if abs(p) < 1e-16 or abs(q) < 1e-16: + # Do not update if probabilities are + # too small + continue + + ratio = exp(log(p) - log(q)) + + proposal.lambdas[child][parent] += eta_rate(index) * ( + ratio - proposal.lambdas[child][parent]) + + return proposal