From 2e42b1f47a1beb6b173eaf5b07c2d878cb1ad48d Mon Sep 17 00:00:00 2001
From: mbaumBielefeld <mbaum@TechFak.Uni-Bielefeld.DE>
Date: Thu, 31 Jan 2013 17:35:55 +0100
Subject: [PATCH] Extended examples/MCMC.py, improved MarkovChainSampler to use
less memory (lukas idea), started working on division
---
examples/MCMC.py | 22 +++++++++++--
primo/reasoning/MCMC.py | 4 +--
primo/reasoning/density/ProbabilityTable.py | 35 +++++++++++++--------
3 files changed, 43 insertions(+), 18 deletions(-)
diff --git a/examples/MCMC.py b/examples/MCMC.py
index 9faa455..840bf49 100644
--- a/examples/MCMC.py
+++ b/examples/MCMC.py
@@ -4,6 +4,7 @@ from primo.core import BayesNet
from primo.reasoning import DiscreteNode
from primo.reasoning import MarkovChainSampler
from primo.reasoning import GibbsTransitionModel
+from primo.reasoning.density import ProbabilityTable
import numpy
#Construct some simple BayesianNetwork
@@ -28,7 +29,22 @@ transition_model = GibbsTransitionModel()
mcs = MarkovChainSampler()
initial_state={burglary:"Safe",alarm:"Silent"}
-chain = mcs.generateMarkovChain(bn, 20, transition_model, initial_state)
+chain = mcs.generateMarkovChain(bn, 1000000, transition_model, initial_state)
+
+#for c in chain:
+# print c
+
+
+pt = ProbabilityTable()
+pt.add_variable(burglary)
+pt.add_variable(alarm)
+pt.to_jpt_by_states(chain)
+print "----joint-probability----"
+print pt
+print "----burglary----"
+print pt.marginalization(alarm)
+print "----alarm----"
+print pt.division(burglary)
+
+bn.draw()
-for c in chain:
- print c
diff --git a/primo/reasoning/MCMC.py b/primo/reasoning/MCMC.py
index 22d8ec5..537a3a6 100644
--- a/primo/reasoning/MCMC.py
+++ b/primo/reasoning/MCMC.py
@@ -1,5 +1,5 @@
import random
-
+import copy
def weighted_random(weights):
counter = random.random() * sum(weights)
for i,w in enumerate(weights):
@@ -57,7 +57,7 @@ class MarkovChainSampler(object):
def generateMarkovChain(self, network, time_steps, transition_model, initial_state):
state=initial_state
- for t in range(time_steps):
+ for t in xrange(time_steps):
yield state
state=transition_model.transition(network, state)
diff --git a/primo/reasoning/density/ProbabilityTable.py b/primo/reasoning/density/ProbabilityTable.py
index c63527d..5aecb53 100644
--- a/primo/reasoning/density/ProbabilityTable.py
+++ b/primo/reasoning/density/ProbabilityTable.py
@@ -10,16 +10,11 @@ class ProbabilityTable(Density):
def __init__(self):
super(ProbabilityTable, self).__init__()
-
- #self.owner = owner
- #self.variables = [owner]
-
- #size_of_range = len(owner.value_range)
- #self.table = numpy.ones(size_of_range) / size_of_range
-
self.variables = []
self.table = numpy.array(0)
+
+
def get_table(self):
return self.table
@@ -43,9 +38,21 @@ class ProbabilityTable(Density):
self.table = table
self.variables = nodes
+ def to_jpt_by_states(self, samples):
+ '''This method uses a list of variable-instantiations to change this nodes internal
+ table to represent a joint probability table constructed from the given samples.
+ The Argument samples is a list of pairs (RandomNode, value).'''
+
+ for state in samples:
+ index = self.get_cpt_index(state.items())
+ self.table[index]=self.table[index]+1
+
+ self.normalize_as_jpt()
+
+
def set_probability(self, value, node_value_pairs):
index = self.get_cpt_index(node_value_pairs)
- self.table[tuple(index)]=value
+ self.table[index]=value
def get_cpt_index(self, node_value_pairs):
nodes, values = zip(*node_value_pairs)
@@ -54,7 +61,7 @@ class ProbabilityTable(Density):
index_in_values_list = nodes.index(node)
value = values[index_in_values_list]
index.append(node.value_range.index(value))
- return index
+ return tuple(index)
def is_normalized_as_cpt(self,owner):
@@ -67,9 +74,9 @@ class ProbabilityTable(Density):
def is_normalized_as_jpt(self):
return numpy.sum(self.table) == 1.0
- def normalized_as_jpt(self):
- '''This method returns a new ProbabilityTable which is a normalized version of this ProbabilityTable (As a joint probability).'''
- return self.table/numpy.sum(self.table)
+ def normalize_as_jpt(self):
+ '''This method normalizes this ProbabilityTable so it represents a valid joint probability table'''
+ self.table=self.table*1.0/numpy.sum(self.table)
def multiplication(self, inputFactor):
'''This method returns a unified ProbabilityTable which contains the variables of both; the inputFactor
@@ -164,7 +171,9 @@ class ProbabilityTable(Density):
def division(self, factor):
- raise Exception("Called unimplemented function")
+ '''Returns a new ProbabilityTable which is the result of dividing this one by the one given
+ with the argument factor'''
+
def __str__(self):
return str(self.table)
--
GitLab