Worked on MCMC examples

examples/MCMC_Continuous.py

@@ -5,15 +5,14 @@ from primo.reasoning.density import LinearBetaParameters
 from primo.reasoning.density import LinearGaussParameters
 from primo.reasoning.density import Gauss
 from primo.reasoning import MCMC
+from primo.reasoning.convergence_test import ConvergenceTestSimpleCounting
 
 from primo.reasoning import EvidenceEqual as EvEqual
-from primo.reasoning import EvidenceIntervall as EvInterval
+from primo.reasoning import EvidenceLower as EvLower
 
 import numpy
 
-#Construct some simple BayesianNetwork. In this example it models
-#the linear relationship between the age of a plant and the height that
-#it has grown up to (+noise)
+#Construct some simple BayesianNetwork. 
 
 #topology
 bn = BayesNet()
@@ -28,17 +27,34 @@ bn.add_edge(age,height)
 bn.add_edge(age,diameter)
 
 #parameterization
-age_parameters=LinearExponentialParameters(1.0,{},2.0)
+#Semantics: Many young plants and the higher the age the lower the probabilty
+#->lambda=2.0
+age_parameters=LinearExponentialParameters(0.0,{})
 age.set_density_parameters(age_parameters)
 
-height_parameters=LinearGaussParameters(0.1,{age:1},0.3)
+#Semantics: plants start at 0.1 meters underground and grow each year by 1 meter.
+#           variance is 0.3 
+height_parameters=LinearGaussParameters(-0.1,{age:1},0.3)
 height.set_density_parameters(height_parameters)
 
-diameter_parameters=LinearBetaParameters(0.01,{age:0.2},2,5)
+#Semantics: At small age: low alpha, high beta -> skew to the left: thin plants
+#            At higher age: high alpha, low beta -> skew to the right: thick plants
+diameter_parameters=LinearBetaParameters(-10.0,{age:4.0},10.0,{age:-4.0})
 diameter.set_density_parameters(diameter_parameters)
 
 
-mcmc_ask=MCMC(bn,1000)
+mcmc_ask=MCMC(bn,100000,convergence_test=ConvergenceTestSimpleCounting(50000))
+
+
+print "------PriorMarginal:------"
+pm=mcmc_ask.calculate_PriorMarginal([age],Gauss)
+print pm
+print "Ground truth: mu=0.5 C=[0.25]"
+#pm=mcmc_ask.calculate_PriorMarginal([height,diameter],Gauss)
+pm=mcmc_ask.calculate_PriorMarginal([height],Gauss)
+print pm
+print ""
+
 
 evidence={age:EvEqual(2)}
 
@@ -47,14 +63,18 @@ evidence={age:EvEqual(2)}
 #poe=mcmc_ask.calculate_PoE(evidence)
 #print poe
 
-print "PosteriorMarginal:"
+print "------PosteriorMarginal:------"
 pm=mcmc_ask.calculate_PosteriorMarginal([age,height],evidence,Gauss)
 #pm=mcmc_ask.calculate_PosteriorMarginal([height],evidence,Gauss)
+print "P(age,height|age=2):"
 print pm
+print "Ground truth: age=2, height=mu:1.9,C=0.3"
+
+print "------PropabilityOfEvidence------"
+poe=mcmc_ask.calculate_PoE({age:EvLower(0.347)})
+print "Probabilty that age is lower than it's median:"
+print "p(age<0.347)="+str(poe)
+print "Ground truth=0.5"
+
+
 
-print "PriorMarginal:"
-pm=mcmc_ask.calculate_PriorMarginal([age],Gauss)
-print pm
-#pm=mcmc_ask.calculate_PriorMarginal([height,diameter],Gauss)
-pm=mcmc_ask.calculate_PriorMarginal([height],Gauss)
-print pm

examples/MCMC_Discrete.py

@@ -11,13 +11,15 @@ import pdb
 #Construct some simple BayesianNetwork
 bn = BayesNet()
 burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
-alarm = DiscreteNode("Alarm", ["Ringing", "Silent","Kaputt"])
+alarm = DiscreteNode("Alarm", ["Ringing", "Silent","Destroyed"])
 
 bn.add_node(burglary)
 bn.add_node(alarm)
 
 bn.add_edge(burglary,alarm)
 
+
+#Parametrize the network
 burglary_cpt=numpy.array([0.2,0.8])
 burglary.set_probability_table(burglary_cpt, [burglary])
 
@@ -25,29 +27,36 @@ alarm_cpt=numpy.array([[0.8,0.15,0.05],[0.05,0.9,0.05]])
 alarm.set_probability_table(alarm_cpt, [burglary,alarm])
 
 
+#Get some inference object
+mcmc_ask=MCMC(bn,5000)
 
-
-mcmc_ask=MCMC(bn,1000)
-
+#Do some Inferences
 evidence={burglary:EvEq("Intruder")}
 
-
-print "---ProbabilityOfEvidence:---" 
+print "-------ProbabilityOfEvidence:-------" 
 poe=mcmc_ask.calculate_PoE(evidence)
-print poe
+print "p(evidence=Intruder)="+str(poe)
+print "Ground truth=0.2\n"
 
-print "---PosteriorMarginal:---"
+print "-------PosteriorMarginal:-------"
 pm=mcmc_ask.calculate_PosteriorMarginal([alarm],evidence,ProbabilityTable)
-print pm
+print "P(alarm|burglary=Intruder)="+str(pm)
+print "Ground truth=[0.8, 0.15, 0.05]\n"
 
-print "---PriorMarginal:---"
+print "-------PriorMarginal:-------"
 pm=mcmc_ask.calculate_PriorMarginal([alarm],ProbabilityTable)
-print "Alarm: " + str(pm)
+print "P(Alarm)= " + str(pm)
+print "Ground truth=[0.2, 0.75, 0.05]\n"
+
 pm=mcmc_ask.calculate_PriorMarginal([burglary],ProbabilityTable)
-print "Burglary: " + str(pm)
+print "P(Burglary)= " + str(pm)
+print "Ground truth=[0.2, 0.8]\n"
 
-print "---MAP:---"
+print "-------MAP:-------"
 hyp=mcmc_ask.calculate_MAP([alarm],evidence,ProbabilityTable)
-print str(evidence) + ": " + str(hyp)
-hyp=mcmc_ask.calculate_MAP([alarm],{},ProbabilityTable)
-print str({}) + ": "+str(hyp)
+print "MAP(alarm|burglary=intruder)=" + str(hyp)
+print "Ground truth=\"Ringing\"\n"
+
+hyp=mcmc_ask.calculate_MAP([burglary,alarm],{},ProbabilityTable)
+print "MAP(burglary,alarm)="+str(hyp)
+print "Ground truth=\"Safe\",\"Silent\"\n"

examples/MCMC_Gauss.py

-from primo.core import BayesNet
-from primo.reasoning import ContinuousNodeFactory
-from primo.reasoning.density import LinearGaussParameters
-from primo.reasoning.density import Gauss
-from primo.reasoning import MCMC
-from primo.reasoning import EvidenceEqual as EvEq
-from primo.reasoning import EvidenceLower as EvLower
-from primo.reasoning import EvidenceIntervall as EvInterval
-
-import numpy
-
-#Construct some simple BayesianNetwork. In this example it models
-#the linear relationship between the age of a plant and the height that
-#it has grown up to (+noise)
-
-#topology
-bn = BayesNet()
-cnf=ContinuousNodeFactory()
-age = cnf.createLinearGaussNode("Plant_age")
-height = cnf.createLinearGaussNode("Plant_height")
-diameter = cnf.createLinearGaussNode("Plant_diameter")
-bn.add_node(age)
-bn.add_node(height)
-bn.add_node(diameter)
-bn.add_edge(age,height)
-bn.add_edge(age,diameter)
-
-
-
-
-#parameterization
-age_b0=4.0
-age_var=2.0
-age_b={}
-age_parameters=LinearGaussParameters(age_b0,age_b,age_var)
-age.set_density_parameters(age_parameters)
-
-height_b0=2.0
-height_var=1.5
-height_b={age:0.5}
-height_parameters=LinearGaussParameters(height_b0,height_b,height_var)
-height.set_density_parameters(height_parameters)
-
-diameter_parameters=LinearGaussParameters(0.1,{age:-0.2},0.1)
-diameter.set_density_parameters(diameter_parameters)
-
-
-mcmc_ask=MCMC(bn,1000)
-
-evidence={age:EvEq(2)}
-
-
-#print "ProbabilityOfEvidence: " 
-#poe=mcmc_ask.calculate_PoE(evidence)
-#print poe
-
-print "PosteriorMarginal:"
-pm=mcmc_ask.calculate_PosteriorMarginal([height,diameter],evidence,Gauss)
-print pm
-
-print "PriorMarginal:"
-pm=mcmc_ask.calculate_PriorMarginal([height,diameter],Gauss)
-print pm
-
-print "ProbabilityOfEvidence:"
-poe=mcmc_ask.calculate_PoE({age:EvLower(4.0)})
-print "p(age<4.0):"+str(poe)+ " GroundTruth:0.5"
-poe=mcmc_ask.calculate_PoE({age:EvInterval(4.0-2.0**0.5,4.0)})
-print "p(4.0-2.0^0.5<age<4.0):"+str(poe)+ " GroundTruth:1/3"
-#print "PriorMarginal:"
-#pm=mcmc_ask.calculate_PriorMarginal([alarm])
-#print "Alarm: " + str(pm)
-#pm=mcmc_ask.calculate_PriorMarginal([burglary])
-#print "Burglary: " + str(pm)

primo/reasoning/ContinuousNode.py

@@ -24,6 +24,9 @@ class ContinuousNode(RandomNode):
     def __str__(self):
         return self.name
         
+    def __repr__(self):
+        return "str(ContinuousNode)"+self.name+")"
+        
     def set_density_parameters(self, density_parameters):
         self.cpd.set_parameters(density_parameters)
         
@@ -50,10 +53,15 @@ class ContinuousNode(RandomNode):
         reimplement it by constructing the integral over the normalvariate in the
         intervalls allowed by the evidence and then generate a sample directly.
         Currently this method has O(inf).'''
-        v=random.normalvariate(x,1.0)
+        std_walk=0.1
+        v=random.normalvariate(x,std_walk)
         if self in evidence.keys():
-            while not evidence[self].is_compatible(v):
-                v=random.normalvariate(x,1.0)
+            while not (evidence[self].is_compatible(v) and self.value_range[0]<v and v<self.value_range[1]):
+                v=random.normalvariate(x,std_walk)
+        else:
+            while not (self.value_range[0]<v and v<self.value_range[1]):
+                v=random.normalvariate(x,std_walk)
+        #print "REEEEEEEEEEEEEEEEEEEEEEEEEEEEEE"+self.name+" " +str(v)
         return v
         
     def sample_global(self, state):

primo/reasoning/DiscreteNode.py

@@ -14,6 +14,12 @@ class DiscreteNode(RandomNode):
         self.value_range = value_range
         self.cpd = ProbabilityTable()
         self.cpd.add_variable(self)
+        
+    def __str__(self):
+        return self.name
+        
+    def __repr__(self):
+        return "DiscreteNode("+self.name+")"
 
     def set_probability(self, value, node_value_pairs):
         self.cpd.set_probability(value, node_value_pairs)

primo/reasoning/density/Gauss.py

@@ -61,9 +61,5 @@ class Gauss(Density):
         return self
         
     def __str__(self):
-        ret= "Gauss:\n"
-        ret=ret+ "--mu--\n"
-        ret=ret+ str(self.mu)+"\n"
-        ret=ret+ "--C--\n"
-        ret=ret+ str(self.C)+"\n"
+        ret= "Gauss(\nmu="+str(self.mu)+"\nC="+str(self.C)+")"
         return ret

primo/reasoning/density/LinearBeta.py

@@ -10,6 +10,7 @@ class LinearBetaParameters():
         self.p0=p0
         self.q0=q0
 
+
 class LinearBeta(Density):
     def __init__(self, node):
         self.p0=1
@@ -18,6 +19,9 @@ class LinearBeta(Density):
         self.q={}
         self.node=node
         
+        self.input_scale=0.1
+        self.output_scale=5.0
+        
     def set_parameters(self,parameters):
         self.p0=parameters.p0
         self.q0=parameters.q0
@@ -38,28 +42,35 @@ class LinearBeta(Density):
     def get_probability(self,value, node_value_pairs):
         #print "Probability to compute for value:"+str(value)
         #Compute the offset for the density and displace the value accordingly
-        p = self.p0
-        q = self.q0
-
-        for node,node_value in node_value_pairs:
-            p = p + self.p[node]*node_value
-            q = q + self.q[node]*node_value
+        #p = self.p0
+        #q = self.q0
 
-        return beta(p, q).pdf(value)
+        #for node,node_value in node_value_pairs:
+        #    p = p + self.p[node]*node_value
+        #    q = q + self.q[node]*node_value
+        #p=1.0/(1.0+math.exp(-p))
+        #q=1.0/(1.0+math.exp(-q))
+        p=self._compute_p_given_parents(dict(node_value_pairs))
+        q=self._compute_q_given_parents(dict(node_value_pairs))
+        #print node_value_pairs
+        #print "beta "+str(p)+" "+str(q)
+        probability = beta(p, q).pdf(value)
+        #print "/beta"
+        return probability
         
     def _compute_p_given_parents(self, state):
         x = self.p0
         for node in self.p.keys():
             if node in state.keys():
                 x = x + self.p[node]*state[node]
-        return 1.0/(1.0+math.exp(-x))
+        return self.output_scale*1.0/(1.0+math.exp(-x*self.input_scale))
         
     def _compute_q_given_parents(self, state):
         x = self.q0
         for node in self.q.keys():
             if node in state.keys():
                 x = x + self.q[node]*state[node]
-        return 1.0/(1.0+math.exp(-x))
+        return self.output_scale*1.0/(1.0+math.exp(-x*self.input_scale))
 
 
     def sample_global(self, state):

primo/reasoning/density/LinearExponential.py

@@ -17,6 +17,9 @@ class LinearExponential(Density):
         self.b0=0
         self.node=node
         
+        self.input_scale=1.0
+        self.output_scale=4.0
+        
     def set_parameters(self,parameters):
         self.b=parameters.b
         self.b0=parameters.b0
@@ -33,10 +36,7 @@ class LinearExponential(Density):
     def get_probability(self,value, node_value_pairs):
         
         #Compute the offset for the density and displace the value accordingly
-        x = self.b0
-        for node,value in node_value_pairs:
-            x = x + self.b[node]*value
-        _lambda=1.0/(1.0+math.exp(-x))
+        _lambda = self._compute_lambda_given_parents(dict(node_value_pairs))
         #print "lambda:"+str(_lambda)
         #Evaluate the displaced density at value
         return _lambda*math.exp(-_lambda*value)
@@ -46,9 +46,11 @@ class LinearExponential(Density):
         for node in self.b.keys():
             if node in state.keys():
                 x = x + self.b[node]*state[node]
-        _lambda=1.0/(1.0+math.exp(-x))
+        _lambda=self.output_scale*1.0/(1.0+math.exp(-x*self.input_scale))
         return _lambda
 
     def sample_global(self,state):
         _lambda=self._compute_lambda_given_parents(state)
-        return random.expovariate(_lambda)
+        sample=random.expovariate(_lambda)
+        #print "EXPO-SAMPLE: "+str(sample)+" at lambda: "+str(_lambda)
+        return sample