Commit c3d5f339 authored by Manuel Baum's avatar Manuel Baum
Browse files

Documentation

parent afc13a69
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.density import ExponentialParameters
from primo.reasoning.density import BetaParameters
from primo.reasoning.density import GaussParameters
from primo.reasoning.density import NDGauss
from primo.reasoning import MCMC
from primo.reasoning.convergence_test import ConvergenceTestSimpleCounting
......@@ -12,34 +12,47 @@ from primo.reasoning import EvidenceLower as EvLower
import numpy
#About this example:
#This example shows how approximate inference can be used query a purely continuous
#bayesian network. At first that network is being constructed and afterwards it
#is passed to an MCMC object that is used to answer several kinds of questions:
#-Prior marginal
#-Posterior marginal
#-Probability of evidence
#-Maximum a-posteriori hypothesis
#Construct some simple BayesianNetwork.
#topology
bn = BayesNet()
cnf=ContinuousNodeFactory()
age = cnf.createLinearExponentialNode("Plant_age")
height = cnf.createLinearGaussNode("Plant_height")
diameter = cnf.createLinearBetaNode("Plant_diameter")
age = cnf.createExponentialNode("Plant_age")
height = cnf.createGaussNode("Plant_height")
diameter = cnf.createBetaNode("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
#Semantics: Many young plants and the higher the age the lower the probabilty
#->lambda=2.0
age_parameters=LinearExponentialParameters(0.0,{})
age_parameters=ExponentialParameters(0.0,{})
age.set_density_parameters(age_parameters)
#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_parameters=GaussParameters(-0.1,{age:1},0.3)
height.set_density_parameters(height_parameters)
#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_parameters=BetaParameters(-10.0,{age:4.0},10.0,{age:-4.0})
diameter.set_density_parameters(diameter_parameters)
......@@ -49,22 +62,22 @@ mcmc_ask=MCMC(bn,1000,convergence_test=ConvergenceTestSimpleCounting(500))
print "------PriorMarginal:------"
pm=mcmc_ask.calculate_PriorMarginal([age],Gauss)
pm=mcmc_ask.calculate_PriorMarginal([age],NDGauss)
print pm
print "Ground truth: mu=0.5 C=[0.25]"
pm=mcmc_ask.calculate_PriorMarginal([height],Gauss)
pm=mcmc_ask.calculate_PriorMarginal([height],NDGauss)
print pm
print ""
print "------PosteriorMarginal:------"
pm=mcmc_ask.calculate_PosteriorMarginal([age,height],{age:EvEqual(2)},Gauss)
pm=mcmc_ask.calculate_PosteriorMarginal([age,height],{age:EvEqual(2)},NDGauss)
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([age,height],{age:EvLower(0.1)},NDGauss)
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"
......@@ -78,7 +91,7 @@ print "Ground truth=0.5"
print ""
print "------MAP------"
map_hypothesis=mcmc_ask.calculate_MAP([height,diameter],{},Gauss)
map_hypothesis=mcmc_ask.calculate_MAP([height,diameter],{},NDGauss)
print map_hypothesis
......
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.density import ExponentialParameters
from primo.reasoning.density import BetaParameters
from primo.reasoning.density import GaussParameters
from primo.reasoning.density import NDGauss
from primo.reasoning import MCMC
from primo.reasoning import EvidenceEqual as EvEqual
......@@ -20,13 +20,13 @@ bn = BayesNet()
cnf=ContinuousNodeFactory()
#create the nodes
age = cnf.createLinearExponentialNode("Age")
sun = cnf.createLinearBetaNode("Sun")
ground= cnf.createLinearGaussNode("Ground")
growth= cnf.createLinearGaussNode("Growth")
height = cnf.createLinearBetaNode("Height")
diameter = cnf.createLinearExponentialNode("Diameter")
children = cnf.createLinearExponentialNode("Children")
age = cnf.createExponentialNode("Age")
sun = cnf.createBetaNode("Sun")
ground= cnf.createGaussNode("Ground")
growth= cnf.createGaussNode("Growth")
height = cnf.createBetaNode("Height")
diameter = cnf.createExponentialNode("Diameter")
children = cnf.createExponentialNode("Children")
#add the nodes to the network
bn.add_node(age)
......@@ -47,19 +47,19 @@ bn.add_edge(height,children)
bn.add_edge(ground,children)
#parameterization
age.set_density_parameters(LinearExponentialParameters(0.1,{}))
age.set_density_parameters(ExponentialParameters(0.1,{}))
sun.set_density_parameters(LinearBetaParameters(2,{},2,{}))
sun.set_density_parameters(BetaParameters(2,{},2,{}))
ground.set_density_parameters(LinearGaussParameters(2.0,{},1.5))
ground.set_density_parameters(GaussParameters(2.0,{},1.5))
growth.set_density_parameters(LinearGaussParameters(0.1,{age:5.0,ground:1.0,sun:4.0},2.5))
growth.set_density_parameters(GaussParameters(0.1,{age:5.0,ground:1.0,sun:4.0},2.5))
height.set_density_parameters(LinearBetaParameters(0.1,{growth:1},0.5,{growth:0.5}))
height.set_density_parameters(BetaParameters(0.1,{growth:1},0.5,{growth:0.5}))
diameter.set_density_parameters(LinearExponentialParameters(0.01,{growth:0.2}))
diameter.set_density_parameters(ExponentialParameters(0.01,{growth:0.2}))
children.set_density_parameters(LinearExponentialParameters(0.1,{ground:1.0,height:1.0}))
children.set_density_parameters(ExponentialParameters(0.1,{ground:1.0,height:1.0}))
mcmc_ask=MCMC(bn,1000)
......@@ -68,13 +68,13 @@ evidence={age:EvEqual(2)}
print "PosteriorMarginal:"
pm=mcmc_ask.calculate_PosteriorMarginal([age,height],evidence,Gauss)
pm=mcmc_ask.calculate_PosteriorMarginal([age,height],evidence,NDGauss)
#pm=mcmc_ask.calculate_PosteriorMarginal([height],evidence,Gauss)
print pm
print "PriorMarginal:"
pm=mcmc_ask.calculate_PriorMarginal([age],Gauss)
pm=mcmc_ask.calculate_PriorMarginal([age],NDGauss)
print pm
#pm=mcmc_ask.calculate_PriorMarginal([height,diameter],Gauss)
pm=mcmc_ask.calculate_PriorMarginal([height],Gauss)
pm=mcmc_ask.calculate_PriorMarginal([height],NDGauss)
print pm
......@@ -5,9 +5,19 @@ from primo.reasoning import DiscreteNode
from primo.reasoning.density import ProbabilityTable
from primo.reasoning import MCMC
from primo.reasoning import EvidenceEqual as EvEq
from primo.reasoning import GibbsTransitionModel
import numpy
import pdb
#About this example:
#This example shows how approximate inference can be used query a purely discrete
#bayesian network. At first that network is being constructed and afterwards it
#is passed to an MCMC object that is used to answer several kinds of questions:
#-Prior marginal
#-Posterior marginal
#-Probability of evidence
#-Maximum a-posteriori hypothesis
#Construct some simple BayesianNetwork
bn = BayesNet()
burglary = DiscreteNode("Burglary", ["Intruder","Safe"])
......@@ -28,7 +38,7 @@ alarm.set_probability_table(alarm_cpt, [burglary,alarm])
#Get some inference object
mcmc_ask=MCMC(bn,5000)
mcmc_ask=MCMC(bn,5000,transition_model=GibbsTransitionModel())
#Do some Inferences
evidence={burglary:EvEq("Intruder")}
......
......@@ -31,14 +31,6 @@ class ContinuousNode(RandomNode):
def set_density_parameters(self, density_parameters):
self.cpd.set_parameters(density_parameters)
# def sample_uniform(self):
# sampled_value = random.uniform(self.value_range[0],self.value_range[1])
# return sampled_value
#
# def sample_proposal(self, x=None):
# return self.cpd.sample_proposal(x)
def sample_local(self, x, evidence):
'''
This method can be used to do a random walk in the domain of this node.
......@@ -71,10 +63,6 @@ class ContinuousNode(RandomNode):
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)
......
from primo.reasoning import ContinuousNode
from primo.reasoning.density import LinearGauss
from primo.reasoning.density import LinearExponential
from primo.reasoning.density import LinearBeta
from primo.reasoning.density import Gauss
from primo.reasoning.density import Exponential
from primo.reasoning.density import Beta
class ContinuousNodeFactory(object):
'''This class offers methods for generating ContinuousNodes'''
def __init__(self):
pass
def createLinearGaussNode(self, name):
def createGaussNode(self, name):
'''
Create a LinearGaussNode with linear dependencies on parents.
@param name: The name of the node.
'''
return self.createContinuousNode(name,(-float("Inf"),float("Inf")),LinearGauss)
return self.createContinuousNode(name,(-float("Inf"),float("Inf")),Gauss)
def createLinearExponentialNode(self, name):
def createExponentialNode(self, name):
'''
Create a LinearExponentialNode with linear dependencies on parents.
@param name: The name of the node.
'''
return self.createContinuousNode(name,(0,float("Inf")),LinearExponential)
return self.createContinuousNode(name,(0,float("Inf")),Exponential)
def createLinearBetaNode(self, name):
def createBetaNode(self, name):
'''
Create a LinearBetaNode with linear dependencies on parents.
@param name: The name of the node.
'''
return self.createContinuousNode(name,(0,1),LinearBeta)
return self.createContinuousNode(name,(0,1),Beta)
def createContinuousNode(self,name,value_range,DensityClass):
'''
......
......@@ -33,7 +33,7 @@ class DiscreteNode(RandomNode):
def is_valid(self):
return self.cpd.is_normalized_as_cpt(self)
def sample_global(self, evidence=None):
def sample_global(self, state, evidence):
if evidence==None or not self in evidence.keys():
compatibles=self.value_range
else:
......@@ -42,7 +42,7 @@ class DiscreteNode(RandomNode):
if evidence[self].is_compatible(v):
compatibles.append(v)
return random.choice(compatibles)
return self.cpd.sample_global(state,self,compatibles)
def sample_local(self, x, evidence=None):
if evidence==None or not self in evidence.keys():
......@@ -53,4 +53,4 @@ class DiscreteNode(RandomNode):
if evidence[self].is_compatible(v):
compatibles.append(v)
return random.choice(compatibles)
return random.choice(compatibles),1.0
......@@ -119,13 +119,6 @@ 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,evidence)
# state[var]=value
# else:
# state[var]=var.sample_global(state)
state[var]=var.sample_global(state, evidence)
return state
......@@ -16,9 +16,18 @@ def weighted_random(weights):
class GibbsTransitionModel(object):
'''
Implements Gibbs-sampling. Can be used to constuct a Markov Chain and is
mainly used by MarkovChainSampler.
mainly used by MarkovChainSampler. This transition model can only be used
if the product of each variable and the variables in it's markov blanket
can be computed in closed form. This is currently only the case for discrete
variables / ProbabilityTables, but could possibly extended to the continuous
setting by assuming gaussian forms for the products or using only classes of
pdfs for which closed forms are computable.
After "Probabilistic Graphical Models, Daphne Koller and Nir Friedman"(p.506)
If the pdf-classes used can not offer this kind of computation you should
use the MetropolisHastingsTransitionModel, as it only requires to compute
a single probability, which can way easier be obtained.
Implemented after "Probabilistic Graphical Models, Daphne Koller and Nir Friedman"(p.506)
'''
def __init__(self):
pass
......@@ -101,6 +110,7 @@ class MetropolisHastingsTransitionModel(object):
for node in nodes_to_resample:
#propose a new value for this variable:
current_value = state[node]
#print 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)
......
......@@ -3,15 +3,54 @@ from primo.reasoning import ContinuousNode
from scipy.stats import beta
import random
import math
class LinearBetaParameters():
class BetaParameters():
'''
This represents the parameters for the beta-density class.
'''
def __init__(self, p0, p, q0, q):
#vector of coefficients for parent values to compute p
self.p=p
#vector of coefficients for parent values to compute q
self.q=q
#bias for p
self.p0=p0
#bias for q
self.q0=q0
class LinearBeta(Density):
class Beta(Density):
'''
This class represents an beta probabilty density. Unfortunately this
is currently a little bulky to use as the parameters for the dependency are
not very transparent. This is how the dependency works:
The parameters p,q for the exponential distribution are computed analogous
as the activation of a perceptron with sigmoid-activation function:
output_scale * sigmoid(input_scale* (b0 + b'state)) where b'state means the dot product between b (a vector
of weights) and state (a vector with a value for each variable that this
density depends on). Here: sigmoid=1/(1+exp(-x))
The parameters output_scale and input_scale can be used to strech or compress
the sigmoid.
The reason for this is that the parameters are required to be >0. And with
linear dependencies on the parents this could in no way be guaranteed.
Why the sigmoid function:
I had to guarantee that the parameters are > 0. As i did not want to
impose any restrictions on the value range of the parents it was necessary
to map the support of the parents values to a valid support for the parameters. In
other (and maybe more correct words): The dependency function to compute
p and q needed to be of the form R^n->]0,inf].
The first function that came to my mind was:
weighted sum of the parents values put into an exponential function. This
caused problems due to the fast growth of the exponential.
For that reason i switched to the sigmoid function that guarantees 0<p,q<1.
And because p,q<1 is not very practical output_scale has been introduced
to scale from ]0,1[ to ]0,output_scale[.
input_scale can be used to strech the sigmoid in input_direction.
'''
def __init__(self, node):
self.p0=1
self.q0=1
......@@ -29,10 +68,6 @@ class LinearBeta(Density):
self.q=parameters.q
def add_variable(self, variable):
'''This method needs some serious reworking: Variables should not be permitted
to be parents because of their value range. Instead it should be evaluated
if they can yield parameters for this distribution that are permitted. This
can in any case happen under bad influence coefficients'''
if( not isinstance(variable,ContinuousNode.ContinuousNode)):
raise Exception("Tried to add Variable as parent, but is not a ContinuousNode")
self.p[variable]=0.0
......@@ -62,6 +97,17 @@ class LinearBeta(Density):
def sample_global(self, state, lower_limit, upper_limit):
'''This method can be used to sample from this distribution. It is necessary that
a value for each parent is specified and it is possible to constrain the
value that is being sampled to some intervall.
@param state: A dict (node->value) that specifies a value for each variable
that this density depends on.
@param lower_limit: The lower limit of the interval that is allowed to be
sampled as value.
@param upper_limit: The upper limit of the interval that is allowed to be
sampled as value.
@returns: The sampled value. A real number.
'''
p=self._compute_p_given_parents(state)
q=self._compute_q_given_parents(state)
......
......@@ -4,20 +4,61 @@ import random
import math
from primo.reasoning import ContinuousNode
class LinearExponentialParameters(object):
class ExponentialParameters(object):
'''
This represents the parameters for the Exponential-density class.
'''
def __init__(self, b0, b):
#lambda a-priori
self.b0=b0
#a dict (node,coefficient) that holds the weights that define the depency on each node
self.b=b
class LinearExponential(Density):
class Exponential(Density):
'''
This class represents an exponential probabilty density. Unfortunately this
is currently a little bulky to use as the parameters for the dependency are
not very transparent. This is how the dependency works:
The parameter lambda for the exponential distribution is computed analogous
as the activation of a perceptron with sigmoid-activation function:
output_scale * sigmoid(input_scale* (b0 + b'state)) where b'state means the dot product between b (a vector
of weights) and state (a vector with a value for each variable that this
density depends on). Here: sigmoid=1/(1+exp(-x))
The parameters output_scale and input_scale can be used to strech or compress
the sigmoid.
The reason for this is that the parameter lambda is required to be >0. And with
linear dependencies on the parents this could in no way be guaranteed.
Why the sigmoid function:
I had to guarantee that the parameter lambda is > 0. As i did not want to
impose any restrictions on the value range of the parents it was necessary
to map the support of the parents values to a valid support for lambda. In
other (and maybe more correct words): The dependency function to compute
lambda needed to be of the form R^n->]0,inf].
The first function that came to my mind was:
weighted sum of the parents values put into an exponential function. This
caused problems due to the fast growth of the exponential.
For that reason i switched to the sigmoid function that guarantees 0<lambda<1.
And because lambda<1 is not very practical output_scale has been introduced
to scale from ]0,1[ to ]0,output_scale[.
input_scale can be used to strech the sigmoid in input_direction.
'''
def __init__(self, node):
#the coefficients for the weighted sum of parent-values
self.b={}
#bias
self.b0=0
#the node that this density is associated with
self.node=node
#scaling coefficient to stretch or compress the sigmoid in input-direction
self.input_scale=1.0
#scaling coefficient to stretch or compress the sigmoid in output-direction
self.output_scale=4.0
def set_parameters(self,parameters):
......@@ -49,6 +90,17 @@ class LinearExponential(Density):
return _lambda
def sample_global(self,state, lower_limit, upper_limit):
'''This method can be used to sample from this distribution. It is necessary that
a value for each parent is specified and it is possible to constrain the
value that is being sampled to some intervall.
@param state: A dict (node->value) that specifies a value for each variable
that this density depends on.
@param lower_limit: The lower limit of the interval that is allowed to be
sampled as value.
@param upper_limit: The upper limit of the interval that is allowed to be
sampled as value.
@returns: The sampled value. A real number.
'''
_lambda=self._compute_lambda_given_parents(state)
distribution=scipy.stats.expon(loc=0,scale=1.0/_lambda)
......
# -*- coding: utf-8 -*-
from primo.reasoning.density import Density
from primo.reasoning.ContinuousNode import ContinuousNode
import numpy
import scipy.stats
import random
class GaussParameters(object):
def __init__(self, mu, cov):
self.mu=mu
self.cov=cov
'''
This represents the parameters for the Gauss-density class.
'''
def __init__(self, b0,b,var):
#offset for the mean of this variable, a priori
self.b0=b0
#a vector that holds one weight for each variable that this density depends on
self.b=b
#the variance
self.var=var
class Gauss(Density):
'''TODO: write doc'''
'''
This class represents a one-dimensional normal distribution with a mean that
depends linear on the parents but with invariant variance.'''
def __init__(self):
def __init__(self,variable):
super(Gauss, self).__init__()
self.mu=numpy.array([0.0])
self.C=numpy.array([[1.0]])
self.variables=[]
self.b0=0#numpy.array([0.0])
self.b={}
self.var=1.0
def set_parameters(self,parameters):
self.set_b0(parameters.b0)
self.set_b(parameters.b)
self.set_var(parameters.var)
def add_variable(self, variable):
v_min,v_max=variable.get_value_range()
if not (v_min>= -float("Inf") and v_max <=float("Inf")):
raise Exception("Tried to add Variable into Gaussian densitiy, but variable had wrong value-range")
self.variables.append(variable)
if not isinstance(variable, ContinuousNode):
raise Exception("Tried to add Variable into Gaussian densitiy, but variable is not continuous")
self.b[variable]=0.0
def get_probability(self, x, node_value_pairs):
reduced_mu = self.b0
for node,value in node_value_pairs:
reduced_mu = reduced_mu + self.b[node]*value
return scipy.stats.norm(reduced_mu, numpy.sqrt(self.var)).pdf(x)
m=len(self.variables)
self.mu.resize([m,1])
self.C.resize((m,m))
self.C[m-1,m-1]=1.0
def set_b(self, variable, b):
if not variable in b.keys():
raise Exception("Tried to set dependency-variable b for a variable that has not yet been added to this variable's dependencies")
self.b[variable]=b
def set_parameters(self,parameters):
self.set_mu(parameters.mu)
self.set_cov(parameters.cov)
def set_b(self, b):
if not set(self.b.keys())==set(b.keys()):
raise Exception("The variables given in the new b do not match the old dependencies of this density")
self.b=b
def set_mu(self, mu):
self.mu=mu
def set_b0(self, b0):
self.b0=b0
def set_cov(self, C):
self.C=C
def set_var(self, var):
self.var=var
def sample(self):
return numpy.random.multivariate_normal(self.mu,self.C)