# Copyright (c) 2012, CyberPoint International, LLC
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
# * Neither the name of the CyberPoint International, LLC nor the
# names of its contributors may be used to endorse or promote products
# derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL CYBERPOINT INTERNATIONAL, LLC BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
'''
This module provides tools to represent and handle Bayesian networks with linear Gaussian conditional probability distributions.
A linear Gaussian distribution means that the node has a continuous range of outcomes, with a normal distribution over those outcomes. This normal distribution can be parameterized by a *mean* and a *variance*. A linear Gaussian means that the variance of the node is fixed, while the mean is a linear function of the outcomes of the node's parents. In math terms, the mean :math:`m(u)` of a node *u* is a linear function of the values :math:`x_1,\\dots,x_n` of its parents, each weighted by some coefficient :math:`\\beta_i`:
.. math::
m(u) = \\beta_0 +\\beta_1 x_1 + \\dots + \\beta_n x_n
Linear Gaussians are simple but widely used in statistical modeling.
'''
import random
import math
import sys
from orderedskeleton import OrderedSkeleton
[docs]class LGBayesianNetwork(OrderedSkeleton):
'''
This class represents a Bayesian network with linear Gaussian CPDs. It contains the attributes *V*, *E*, and *Vdata*, as well as the method *randomsample*.
'''
def __init__(self, orderedskeleton=None, nodedata=None):
'''
This class can be called either with or without arguments. If it is called without arguments, none of its attributes are instantiated and it is left to the user to instantiate them manually. If it is called with arguments, the attributes will be loaded directly from the inputs. The arguments must be (in order):
1. *orderedskeleton* -- An instance of the :doc:`OrderedSkeleton <orderedskeleton>` or :doc:`GraphSkeleton <graphskeleton>` (as long as it's ordered) class.
2. *nodedata* -- An instance of the :doc:`NodeData <nodedata>` class.
If these arguments are present, all attributes of the class (*V*, *E*, and *Vdata*) will be automatically copied from the graph skeleton and node data inputs.
This class requires that the *Vdata* attribute gets loaded with a dictionary with node data of the following fomat::
"vertex": {
"parents": ["<name of parent 1>", ... , "<name of parent n>"],
"children": ["<name of child 1>", ... , "<name of child n>"],
"mean_base": <the base mean of the Gaussian distribution>,
"mean_scal": [<scalar for parent 1 outcome>, ... , <scalar for parent n outcome>],
"variance": <variance of the Gaussian distibution>
}
Note that additional keys are possible in the dict of each vertex.
Upon loading, the class will also check that the keys of *Vdata* correspond to the vertices in *V*.
'''
if (orderedskeleton != None and nodedata != None):
try:
self.V = orderedskeleton.V
'''A list of the names of the vertices.'''
self.E = orderedskeleton.E
'''A list of [origin, destination] pairs of vertices that make edges.'''
self.Vdata = nodedata.Vdata
'''A dictionary containing CPD data for the nodes.'''
except:
raise Exception, "Inputs were malformed; first arg must contain V and E attributes and second arg must contain Vdata attribute."
assert sorted(self.V) == sorted(self.Vdata.keys()), "Vertices did not match node data"
[docs] def randomsample(self, n, evidence=None):
'''
Produce *n* random samples from the Bayesian Network and return them in a list.
See above for how the means of linear Gaussians are calculated during sampling.
This function takes the following arguments:
1. *n* -- The number of random samples to produce.
2. *evidence* -- (Optional) A dict containing (vertex: value) pairs that describe the evidence. To be used carefully because it does manually overrides the nodes with evidence instead of affecting the joint probability distribution of the entire graph.
And returns:
A list of *n* independent random samples, each element of which is a dict containing (vertex: value) pairs.
Usage example: this would generate a sequence of 10 random samples::
import json
from libpgm.nodedata import NodeData
from libpgm.graphskeleton import GraphSkeleton
from libpgm.lgbayesiannetwork import LGBayesianNetwork
# load nodedata and graphskeleton
nd = NodeData()
skel = GraphSkeleton()
nd.load("../tests/unittestlgdict.txt") # an input file
skel.load("../tests/unittestdict.txt")
# topologically order graphskeleton
skel.toporder()
# load bayesian network
lgbn = LGBayesianNetwork(skel, nd)
# sample
result = lgbn.randomsample(10)
# output
print json.dumps(result, indent=2)
'''
assert (isinstance(n, int) and n > 0), "Argument must be a positive integer."
random.seed()
seq = []
for _ in range(n):
outcome = dict()
for vertex in self.V:
outcome[vertex] = "default"
def assignnode(s):
if (evidence != None):
if s in evidence.keys():
return evidence[s]
# calculate Bayesian parameters (mean and variance)
mean = self.Vdata[s]["mean_base"]
if (self.Vdata[s]["parents"] != None):
for x in range(len(self.Vdata[s]["parents"])):
parent = self.Vdata[s]["parents"][x]
assert outcome[parent] != 'default', "Graph skeleton was not topologically ordered."
mean += outcome[parent] * self.Vdata[s]["mean_scal"][x]
variance = self.Vdata[s]["variance"]
# draw random outcome from Gaussian (I love python)
return random.gauss(mean, math.sqrt(variance))
for s in self.V:
if (outcome[s] == "default"):
outcome[s] = assignnode(s)
seq.append(outcome)
return seq