Commit 282527ea authored by Manuel Baum's avatar Manuel Baum
Browse files

examples/MCMC_Continuous.py läuft mit exponentialverteilung->linearGauss,...

examples/MCMC_Continuous.py läuft mit exponentialverteilung->linearGauss, Evidenz als Klassen und Kontinuierliche PoE läuft
parent 3bacac8d
......@@ -4,6 +4,7 @@ from primo.core import BayesNet
from primo.reasoning import DiscreteNode
from primo.reasoning.density import ProbabilityTable
from primo.reasoning import MCMC
from primo.reasoning import EvidenceEqual as EvEq
import numpy
import pdb
......@@ -28,7 +29,7 @@ alarm.set_probability_table(alarm_cpt, [burglary,alarm])
mcmc_ask=MCMC(bn)
evidence={burglary:"Intruder"}
evidence={burglary:EvEq("Intruder")}
print "---ProbabilityOfEvidence:---"
......
......@@ -2,9 +2,13 @@ from primo.core import BayesNet
from primo.reasoning import ContinuousNodeFactory
from primo.reasoning.density import LinearExponentialParameters
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 import EvidenceEqual as EvEqual
from primo.reasoning import EvidenceIntervall as EvInterval
import numpy
#Construct some simple BayesianNetwork. In this example it models
......@@ -15,7 +19,7 @@ import numpy
bn = BayesNet()
cnf=ContinuousNodeFactory()
age = cnf.createLinearExponentialNode("Plant_age")
height = cnf.createLinearExponentialNode("Plant_height")
height = cnf.createLinearGaussNode("Plant_height")
#diameter = cnf.createLinearBetaNode("Plant_diameter")
bn.add_node(age)
bn.add_node(height)
......@@ -27,7 +31,7 @@ bn.add_edge(age,height)
age_parameters=LinearExponentialParameters(4,{})
age.set_density_parameters(age_parameters)
height_parameters=LinearExponentialParameters(0.1,{age:1})
height_parameters=LinearGaussParameters(0.1,{age:1},0.3)
height.set_density_parameters(height_parameters)
#diameter_parameters=LinearGaussParameters(0.1,{age:-0.2},0.1)
......@@ -36,7 +40,7 @@ height.set_density_parameters(height_parameters)
mcmc_ask=MCMC(bn)
evidence={age:2}
evidence={age:EvEqual(2)}
#print "ProbabilityOfEvidence: "
......@@ -44,18 +48,13 @@ evidence={age:2}
#print poe
print "PosteriorMarginal:"
#pm=mcmc_ask.calculate_PosteriorMarginal([height,diameter],evidence,Gauss)
pm=mcmc_ask.calculate_PosteriorMarginal([age,height],evidence,Gauss)
#pm=mcmc_ask.calculate_PosteriorMarginal([height],evidence,Gauss)
#print pm
print pm
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)
pm=mcmc_ask.calculate_PriorMarginal([height],Gauss)
print pm
#print "PriorMarginal:"
#pm=mcmc_ask.calculate_PriorMarginal([alarm])
#print "Alarm: " + str(pm)
#pm=mcmc_ask.calculate_PriorMarginal([burglary])
#print "Burglary: " + str(pm)
......@@ -3,6 +3,9 @@ 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
......@@ -44,7 +47,7 @@ diameter.set_density_parameters(diameter_parameters)
mcmc_ask=MCMC(bn)
evidence={age:2}
evidence={age:EvEq(2)}
#print "ProbabilityOfEvidence: "
......@@ -59,6 +62,11 @@ 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)
......
......@@ -41,3 +41,9 @@ class Node(object):
def __str__(self):
print self.name
return self.name
def sample_gobal(self, evidence=None):
raise Exception("Called unimplemented Method")
def sample_local(self, x, evidence=None):
raise Exception("Called unimplemented Method")
......@@ -19,6 +19,29 @@ class ContinuousNode(Node):
def sample_proposal(self, x=None):
return self.cpd.sample_proposal(x)
def sample_local(self, x, evidence=None):
'''This is the most simple and stupid implementation of the method. It
uses bogo-search to find a sample that fits the evidence. You could
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)
if evidence != None:
while not evidence.is_compatible(v):
v=random.normalvariate(x,1.0)
return v
def sample_global(self, evidence=None):
'''Simple, Stupid and O(inf). Improvement idea see comment on sample_local()'''
if evidence==None:
return self.cpd.sample_global()
else:
proposal=self.cpd.sample_global()
while not evidence.is_compatible(proposal):
proposal=self.cpd.sample_global()
return proposal
def get_probability(self, value, node_value_pairs):
return self.cpd.get_probability(value, node_value_pairs)
......@@ -27,8 +27,24 @@ class DiscreteNode(RandomNode):
def is_valid(self):
return self.cpd.is_normalized_as_cpt(self)
def sample_uniform(self):
return random.choice(self.value_range)
def sample_global(self, evidence=None):
if evidence==None:
compatibles=self.value_range
else:
compatibles=[]
for v in self.value_range:
if evidence.is_compatible(v):
compatibles.append(v)
def sample_proposal(self, current_value=None):
return random.choice(self.value_range)
return random.choice(compatibles)
def sample_local(self, x, evidence=None):
if evidence==None:
compatibles=self.value_range
else:
compatibles=[]
for v in self.value_range:
if evidence.is_compatible(v):
compatibles.append(v)
return random.choice(compatibles)
class Evidence(object):
def __init__(self):
pass
def is_compatible(self, value):
raise Exception("Not defined for this kind of Evidence")
def get_unambigous_value(self):
return None
class EvidenceEqual(Evidence):
def __init__(self, value):
self.value=value
def is_compatible(self, value):
return self.value==value
def get_unambigous_value(self):
return self.value
class EvidenceIntervall(Evidence):
def __init__(self,min_val,max_val):
self.min_val=min_val
self.max_val=max_val
def is_compatible(self, value):
return self.min_val <= value and value<=self.max_val
class EvidenceLower(Evidence):
def __init__(self,limit):
self.limit=limit
def is_compatible(self, value):
return value<self.limit
class EvidenceHigher(Evidence):
def __init__(self,limit):
self.limit=limit
def is_compatible(self, value):
return value>self.limit
......@@ -59,9 +59,9 @@ class MCMC(object):
number_of_samples=0
for state in chain:
compatible = True
for node,value in evidence.items():
for node,node_evidence in evidence.items():
if state[node] != value:
if not node_evidence.is_compatible(state[node]):
compatible = False
break
......@@ -77,8 +77,12 @@ class MCMC(object):
state=[]
for var in self.bn.get_nodes([]):
if var in evidence.keys():
state.append((var,evidence[var]))
unambigous=evidence[var].get_unambigous_value()
if unambigous == None:
state.append((var,var.sample_global(evidence[var])))
else:
state.append((var,unambigous))
else:
state.append((var,var.sample_proposal()))
state.append((var,var.sample_global()))
return dict(state)
......@@ -75,7 +75,7 @@ class MetropolisHastingsTransitionModel(object):
for node in nodes_to_resample:
#propose a new value for this variable:
current_value = state[node]
proposed_value = node.sample_proposal(current_value)
proposed_value = node.sample_local(current_value)
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)
......@@ -102,7 +102,7 @@ class MarkovChainSampler(object):
state=initial_state
if evidence:
for node in evidence.keys():
if state[node] != evidence[node]:
if not evidence[node].is_compatible(state[node]):
raise Exception("The evidence given does not fit to the initial_state specified")
constant_nodes = evidence.keys()
else:
......@@ -113,6 +113,7 @@ class MarkovChainSampler(object):
state=transition_model.transition(network, state, constant_nodes)
#finally sample from the target distribution
for t in xrange(time_steps):
#print state
if variables_of_interest:
yield self._reduce_state_to_variables_of_interest(state, variables_of_interest)
else:
......
......@@ -3,6 +3,8 @@ from DiscreteNode import DiscreteNode
from ContinuousNode import ContinuousNode
from ContinuousNodeFactory import ContinuousNodeFactory
from Evidence import EvidenceEqual, EvidenceIntervall, EvidenceLower, EvidenceHigher
from MarkovChainSampler import MarkovChainSampler
from MarkovChainSampler import GibbsTransitionModel
from MarkovChainSampler import MetropolisHastingsTransitionModel
......
......@@ -5,8 +5,9 @@ import math
class LinearExponentialParameters(object):
def __init__(self, b0, b):
self.b=b
self.b0=b0
self.b=b
class LinearExponential(Density):
......@@ -29,13 +30,15 @@ class LinearExponential(Density):
self.b[variable]=0.0
def get_probability(self,value, node_value_pairs):
if value<=0:
return 0
lambda_parameter = self.b0
for node,value in node_value_pairs:
lambda_parameter = lambda_parameter + self.b[node]*value
#print scipy.stats.expon(lambda_parameter).pdf(value)
return lambda_parameter*math.exp(-lambda_parameter*value)
def sample_proposal(self, x=None):
if x==None:
x=self.b0
return random.expovariate(x)
def sample_global(self):
print self.b0
return random.expovariate(self.b0)
# -*- coding: utf-8 -*-
from primo.reasoning.density import Density
from primo.reasoning.ContinuousNode import ContinuousNode
import numpy
import scipy.stats
import random
class LinearGaussParameters(object):
def __init__(self, b0,b,var):
......@@ -26,8 +28,9 @@ class LinearGauss(Density):
self.set_var(parameters.var)
def add_variable(self, variable):
if( not variable.get_value_range() == (-float("Inf"),float("Inf"))):
raise Exception("Tried to add Variable into Gaussian densitiy, but variable had wrong value-range")
if not isinstance(variable, ContinuousNode):
raise Exception("Tried to add Variable into Gaussian densitiy, but variable is not continuous")
self.b[variable]=0.0
......@@ -55,8 +58,6 @@ class LinearGauss(Density):
def set_var(self, var):
self.var=var
def sample_proposal(self, x=None):
if x==None:
x=self.b0
return random.normalvariate(x,1.0)
def sample_global(self):
return random.normalvariate(self.b0,self.var**0.5)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment