It runs now but results are still weird

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):
     """

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