fixed some bug in computation of acceptance probability in MetropolisHastings,...

fixed some bug in computation of acceptance probability in MetropolisHastings, now works with restrained value ranges too. Also changed sample_local in continuous node to direct sampling so there is no rejection sampling anymore

examples/MCMC_Continuous.py

@@ -43,10 +43,12 @@ diameter_parameters=LinearBetaParameters(-10.0,{age:4.0},10.0,{age:-4.0})
 diameter.set_density_parameters(diameter_parameters)
 
 
-mcmc_ask=MCMC(bn,100000,convergence_test=ConvergenceTestSimpleCounting(50000))
+mcmc_ask=MCMC(bn,10000,convergence_test=ConvergenceTestSimpleCounting(5000))
 
 
 print "------PriorMarginal:------"
+
+
 pm=mcmc_ask.calculate_PriorMarginal([age],Gauss)
 print pm
 print "Ground truth: mu=0.5 C=[0.25]"
@@ -56,19 +58,20 @@ print pm
 print ""
 
 
-evidence={age:EvEqual(2)}
-
-
-#print "ProbabilityOfEvidence: " 
-#poe=mcmc_ask.calculate_PoE(evidence)
-#print poe
-
 print "------PosteriorMarginal:------"
-pm=mcmc_ask.calculate_PosteriorMarginal([age,height],evidence,Gauss)
+pm=mcmc_ask.calculate_PosteriorMarginal([age,height],{age:EvEqual(2)},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 ""
+
+pm=mcmc_ask.calculate_PosteriorMarginal([age,height],{age:EvLower(0.1)},Gauss)
+#pm=mcmc_ask.calculate_PosteriorMarginal([height],evidence,Gauss)
+print "P(age,height|age<0.1):"
+print pm
+print "Ground truth: age=0:0.1, height=mu:-0.1:0.0,C=0.3"
+print ""
 
 print "------PropabilityOfEvidence------"
 poe=mcmc_ask.calculate_PoE({age:EvLower(0.347)})

primo/reasoning/ContinuousNode.py

@@ -2,6 +2,7 @@
 
 from primo.reasoning import RandomNode
 import random
+import scipy
 
 class ContinuousNode(RandomNode):
     '''
@@ -53,18 +54,44 @@ 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).'''
-        std_walk=0.1
-        v=random.normalvariate(x,std_walk)
+        std_walk=1.0
+        
+        #intersect possible evidence-interval with value_range:
         if self in evidence.keys():
-            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)
+            evidence_range=evidence[self].get_interval()
+            lower_limit=max(self.value_range[0],evidence_range[0])
+            upper_limit=min(self.value_range[1],evidence_range[1])
         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
+            lower_limit=self.value_range[0]
+            upper_limit=self.value_range[1]
+        
+        if lower_limit==upper_limit:
+            v=lower_limit
+        if lower_limit>upper_limit:
+            raise Exception("Intersection of random variable's value_range and"
+                "allowed Interval for Evidence is empty - no sampling possible")
+        
+        #v=random.normalvariate(x,std_walk)
+        #while not (lower_limit<=v and v<=upper_limit):
+        #    v=random.normalvariate(x,std_walk)
+        
+        #generate the actual sample
+        distribution=scipy.stats.norm(x, std_walk)
+        lower_cdf=distribution.cdf(lower_limit)
+        upper_cdf=distribution.cdf(upper_limit)
+        
+        sample_in_integral=random.uniform(lower_cdf, upper_cdf)
+        
+        sample=distribution.ppf(sample_in_integral)
         
-    def sample_global(self, state):
+        
+        
+        a=scipy.stats.norm(self.value_range[0], std_walk).cdf(x)
+        b=scipy.stats.norm(self.value_range[0], std_walk).cdf(sample)
+        cdf_ratio = a/b 
+        return sample,cdf_ratio
+        
+    def sample_global(self, state, evidence):
         '''
         This method can be used to sample from this local distribution.
         
@@ -72,7 +99,31 @@ class ContinuousNode(RandomNode):
             values of this nodes parents in this dict and the conditional 
             probability density will be adjusted accordingly.
         '''
-        proposal=self.cpd.sample_global(state)
+        #is there some evidence for this node?
+        if self in evidence.keys():
+            #if only one value is allowed we can return it immediatly
+            unique=evidence[self].get_unique_value()
+            if unique!=None:
+                return unique
+            #if a whole interval is allowed intersect it with this variable's
+            #value_range to get limits for sampling
+            else:
+                evidence_range=evidence[self].get_interval()
+                lower_limit=max(self.value_range[0],evidence_range[0])
+                upper_limit=min(self.value_range[1],evidence_range[1])
+        #without evidence this variable's value_range represents limits for sampling
+        else:
+            lower_limit=self.value_range[0]
+            upper_limit=self.value_range[1]
+        #check if only one value is allowed and in case return immediatly
+        if lower_limit==upper_limit:
+            return lower_limit
+        #check for empty interval
+        if lower_limit>upper_limit:
+            raise Exception("Intersection of random variable's value_range and"
+                "allowed Interval for Evidence is empty - no sampling possible")
+        
+        proposal=self.cpd.sample_global(state,lower_limit,upper_limit)
         return proposal
         
     def get_probability(self, value, state):

primo/reasoning/Evidence.py

@@ -41,7 +41,7 @@ class EvidenceEqual(Evidence):
     def get_unique_value(self):
         return self.value
         
-class EvidenceIntervall(Evidence):
+class EvidenceInterval(Evidence):
     '''
     This class can be used to specify evidence that a variable has taken on
     some value in a defined interval.
@@ -54,27 +54,26 @@ class EvidenceIntervall(Evidence):
     def is_compatible(self, value):
         return self.min_val <= value and value<=self.max_val
         
+    def get_interval(self):
+        return self.min_val,self.max_val
+        
 
-class EvidenceLower(Evidence):
+class EvidenceLower(EvidenceInterval):
     '''
     This class can be used to specify evidence that a variable has taken on
     some value lower than some threshold.
     e.g. a<3
     '''
     def __init__(self,limit):
-        self.limit=limit
-        
-    def is_compatible(self, value):
-        return value<self.limit
+        super(EvidenceLower, self).__init__(float("-inf"),limit)
+
         
-class EvidenceHigher(Evidence):
+class EvidenceHigher(EvidenceInterval):
     '''
     This class can be used to specify evidence that a variable has taken on
     some value higher than some threshold.
     e.g. a>3
     '''
     def __init__(self,limit):
-        self.limit=limit
-        
-    def is_compatible(self, value):
-        return value>self.limit
+        super(EvidenceLower, self).__init__(limit,float("inf"))
+

primo/reasoning/MCMC.py

@@ -119,12 +119,13 @@ class MCMC(object):
         '''
         state={}
         for var in self.bn.get_nodes_in_topological_sort():
-            if var in evidence.keys():
-                value=evidence[var].get_unique_value()
-                if value == None:
-                    value=var.sample_global(state)
-                state[var]=value
-            else:
-                state[var]=var.sample_global(state)
+#            if var in evidence.keys():
+#                value=evidence[var].get_unique_value()
+#                if value == None:
+#                    value=var.sample_global(state,evidence)
+#                state[var]=value
+#            else:
+#                state[var]=var.sample_global(state)
+            state[var]=var.sample_global(state, evidence)
         return state
 

primo/reasoning/MarkovChainSampler.py

@@ -101,12 +101,12 @@ class MetropolisHastingsTransitionModel(object):
         for node in nodes_to_resample:
             #propose a new value for this variable:
             current_value = state[node]
-            proposed_value = node.sample_local(current_value, extern_evidence)
+            proposed_value, cdf_ratio = node.sample_local(current_value, extern_evidence)
             
             p_of_proposal_given_mb = self._compute_p_of_value_given_mb(network, state, node, proposed_value)
             p_of_current_given_mb = self._compute_p_of_value_given_mb(network, state, node, current_value)
-
-            acceptance_probability = min(1.0,  p_of_proposal_given_mb/p_of_current_given_mb * 1.0/1.0)
+            #print "acceptance_probability = min(1.0,  "+str(p_of_proposal_given_mb)+" / "+str(p_of_current_given_mb) + " * "+str(cdf_ratio)
+            acceptance_probability = min(1.0,  p_of_proposal_given_mb/p_of_current_given_mb * cdf_ratio * 1.0/1.0)
             if random.random() <= acceptance_probability:
                 state[node]=proposed_value
             

primo/reasoning/__init__.py

@@ -3,7 +3,7 @@ from DiscreteNode import DiscreteNode
 from ContinuousNode import ContinuousNode
 from ContinuousNodeFactory import ContinuousNodeFactory
 
-from Evidence import EvidenceEqual, EvidenceIntervall, EvidenceLower, EvidenceHigher
+from Evidence import EvidenceEqual, EvidenceInterval, EvidenceLower, EvidenceHigher
 
 from MarkovChainSampler import MarkovChainSampler
 from MarkovChainSampler import GibbsTransitionModel

primo/reasoning/density/LinearBeta.py

@@ -73,9 +73,31 @@ class LinearBeta(Density):
         return self.output_scale*1.0/(1.0+math.exp(-x*self.input_scale))
 
 
-    def sample_global(self, state):
+    def sample_global(self, state, lower_limit, upper_limit):
         p=self._compute_p_given_parents(state)
         q=self._compute_q_given_parents(state)
-        value=random.betavariate(p,q)
+        #value=random.betavariate(p,q)
         #print "Sampled:"+str(value)
-        return value
+        #return value
+        
+        
+        
+        
+        
+        
+        #_lambda=self._compute_lambda_given_parents(state)
+        distribution=beta(p,q)
+        
+        lower_cdf=distribution.cdf(lower_limit)
+        upper_cdf=distribution.cdf(upper_limit)
+        
+        sample_in_integral=random.uniform(lower_cdf, upper_cdf)
+        
+        sample=distribution.ppf(sample_in_integral)
+        return sample
+        
+        
+        
+        
+        
+        

primo/reasoning/density/LinearExponential.py

@@ -37,7 +37,6 @@ class LinearExponential(Density):
         
         #Compute the offset for the density and displace the value accordingly
         _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)
 
@@ -49,8 +48,17 @@ class LinearExponential(Density):
         _lambda=self.output_scale*1.0/(1.0+math.exp(-x*self.input_scale))
         return _lambda
 
-    def sample_global(self,state):
+    def sample_global(self,state, lower_limit, upper_limit):
         _lambda=self._compute_lambda_given_parents(state)
-        sample=random.expovariate(_lambda)
+        distribution=scipy.stats.expon(loc=0,scale=1.0/_lambda)
+        
+        lower_cdf=distribution.cdf(lower_limit)
+        upper_cdf=distribution.cdf(upper_limit)
+        
+        sample_in_integral=random.uniform(lower_cdf, upper_cdf)
+        
+        sample=distribution.ppf(sample_in_integral)
+        
+        #sample=random.expovariate(_lambda)
         #print "EXPO-SAMPLE: "+str(sample)+" at lambda: "+str(_lambda)
         return sample

primo/reasoning/density/LinearGauss.py

@@ -65,6 +65,17 @@ class LinearGauss(Density):
                 x = x + self.b[node]*state[node]
         return x
         
-    def sample_global(self,state):
-        return random.normalvariate(self._compute_offset_given_parents(state),self.var**0.5)
+    def sample_global(self,state,lower_limit,upper_limit):
+        
+        distribution=scipy.stats.norm(self._compute_offset_given_parents(state), self.var**0.5)
+        
+        lower_cdf=distribution.cdf(lower_limit)
+        upper_cdf=distribution.cdf(upper_limit)
+        
+        sample_in_integral=random.uniform(lower_cdf, upper_cdf)
+        
+        sample=distribution.ppf(sample_in_integral)
+    
+    
+        return sample