From 3ca3d5905aa179fcf8e17fb68f9360b26f0f9f53 Mon Sep 17 00:00:00 2001 From: Manuel Baum <mbaum@techfak.uni-bielefeld.de> Date: Tue, 9 Jul 2013 17:09:07 +0200 Subject: [PATCH] Worked on MCMC examples --- examples/MCMC_Continuous.py | 50 +++++++++---- examples/MCMC_Discrete.py | 41 ++++++----- examples/MCMC_Gauss.py | 74 -------------------- primo/reasoning/ContinuousNode.py | 14 +++- primo/reasoning/DiscreteNode.py | 6 ++ primo/reasoning/density/Gauss.py | 6 +- primo/reasoning/density/LinearBeta.py | 29 +++++--- primo/reasoning/density/LinearExponential.py | 14 ++-- 8 files changed, 106 insertions(+), 128 deletions(-) delete mode 100644 examples/MCMC_Gauss.py diff --git a/examples/MCMC_Continuous.py b/examples/MCMC_Continuous.py index 4ca269b..b419255 100644 --- a/examples/MCMC_Continuous.py +++ b/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 diff --git a/examples/MCMC_Discrete.py b/examples/MCMC_Discrete.py index 0a98ec2..010710a 100644 --- a/examples/MCMC_Discrete.py +++ b/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" diff --git a/examples/MCMC_Gauss.py b/examples/MCMC_Gauss.py deleted file mode 100644 index 21e59ff..0000000 --- a/examples/MCMC_Gauss.py +++ /dev/null @@ -1,74 +0,0 @@ -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) diff --git a/primo/reasoning/ContinuousNode.py b/primo/reasoning/ContinuousNode.py index 77149f3..10c6ce1 100644 --- a/primo/reasoning/ContinuousNode.py +++ b/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): diff --git a/primo/reasoning/DiscreteNode.py b/primo/reasoning/DiscreteNode.py index 90bf85f..551132e 100644 --- a/primo/reasoning/DiscreteNode.py +++ b/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) diff --git a/primo/reasoning/density/Gauss.py b/primo/reasoning/density/Gauss.py index 6e9c90c..3bb5587 100644 --- a/primo/reasoning/density/Gauss.py +++ b/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 diff --git a/primo/reasoning/density/LinearBeta.py b/primo/reasoning/density/LinearBeta.py index f141dea..d1bc190 100644 --- a/primo/reasoning/density/LinearBeta.py +++ b/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): diff --git a/primo/reasoning/density/LinearExponential.py b/primo/reasoning/density/LinearExponential.py index 2e6e464..54c2e18 100644 --- a/primo/reasoning/density/LinearExponential.py +++ b/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 -- GitLab