"""Probability models (Chapter 13-15)"""
from collections import defaultdict
from functools import reduce
from agents import Agent
from utils import *
def DTAgentProgram(belief_state):
"""
[Figure 13.1]
A decision-theoretic agent.
"""
def program(percept):
belief_state.observe(program.action, percept)
program.action = max(belief_state.actions(), key=belief_state.expected_outcome_utility)
return program.action
program.action = None
return program
# ______________________________________________________________________________
class ProbDist:
"""A discrete probability distribution. You name the random variable
in the constructor, then assign and query probability of values.
>>> P = ProbDist('Flip'); P['H'], P['T'] = 0.25, 0.75; P['H']
0.25
>>> P = ProbDist('X', {'lo': 125, 'med': 375, 'hi': 500})
>>> P['lo'], P['med'], P['hi']
(0.125, 0.375, 0.5)
"""
def __init__(self, var_name='?', freq=None):
"""If freq is given, it is a dictionary of values - frequency pairs,
then ProbDist is normalized."""
self.prob = {}
self.var_name = var_name
self.values = []
if freq:
for (v, p) in freq.items():
self[v] = p
self.normalize()
def __getitem__(self, val):
"""Given a value, return P(value)."""
try:
return self.prob[val]
except KeyError:
return 0
def __setitem__(self, val, p):
"""Set P(val) = p."""
if val not in self.values:
self.values.append(val)
self.prob[val] = p
def normalize(self):
"""Make sure the probabilities of all values sum to 1.
Returns the normalized distribution.
Raises a ZeroDivisionError if the sum of the values is 0."""
total = sum(self.prob.values())
if not np.isclose(total, 1.0):
for val in self.prob:
self.prob[val] /= total
return self
def show_approx(self, numfmt='{:.3g}'):
"""Show the probabilities rounded and sorted by key, for the
sake of portable doctests."""
return ', '.join([('{}: ' + numfmt).format(v, p) for (v, p) in sorted(self.prob.items())])
def __repr__(self):
return "P({})".format(self.var_name)
class JointProbDist(ProbDist):
"""A discrete probability distribute over a set of variables.
>>> P = JointProbDist(['X', 'Y']); P[1, 1] = 0.25
>>> P[1, 1]
0.25
>>> P[dict(X=0, Y=1)] = 0.5
>>> P[dict(X=0, Y=1)]
0.5"""
def __init__(self, variables):
self.prob = {}
self.variables = variables
self.vals = defaultdict(list)
def __getitem__(self, values):
"""Given a tuple or dict of values, return P(values)."""
values = event_values(values, self.variables)
return ProbDist.__getitem__(self, values)
def __setitem__(self, values, p):
"""Set P(values) = p. Values can be a tuple or a dict; it must
have a value for each of the variables in the joint. Also keep track
of the values we have seen so far for each variable."""
values = event_values(values, self.variables)
self.prob[values] = p
for var, val in zip(self.variables, values):
if val not in self.vals[var]:
self.vals[var].append(val)
def values(self, var):
"""Return the set of possible values for a variable."""
return self.vals[var]
def __repr__(self):
return "P({})".format(self.variables)
def event_values(event, variables):
"""Return a tuple of the values of variables in event.
>>> event_values ({'A': 10, 'B': 9, 'C': 8}, ['C', 'A'])
(8, 10)
>>> event_values ((1, 2), ['C', 'A'])
(1, 2)
"""
if isinstance(event, tuple) and len(event) == len(variables):
return event
else:
return tuple([event[var] for var in variables])
# ______________________________________________________________________________
def enumerate_joint_ask(X, e, P):
"""
[Section 13.3]
Return a probability distribution over the values of the variable X,
given the {var:val} observations e, in the JointProbDist P.
>>> P = JointProbDist(['X', 'Y'])
>>> P[0,0] = 0.25; P[0,1] = 0.5; P[1,1] = P[2,1] = 0.125
>>> enumerate_joint_ask('X', dict(Y=1), P).show_approx()
'0: 0.667, 1: 0.167, 2: 0.167'
"""
assert X not in e, "Query variable must be distinct from evidence"
Q = ProbDist(X) # probability distribution for X, initially empty
Y = [v for v in P.variables if v != X and v not in e] # hidden variables.
for xi in P.values(X):
Q[xi] = enumerate_joint(Y, extend(e, X, xi), P)
return Q.normalize()
def enumerate_joint(variables, e, P):
"""Return the sum of those entries in P consistent with e,
provided variables is P's remaining variables (the ones not in e)."""
if not variables:
return P[e]
Y, rest = variables[0], variables[1:]
return sum([enumerate_joint(rest, extend(e, Y, y), P) for y in P.values(Y)])
# ______________________________________________________________________________
class BayesNet:
"""Bayesian network containing only boolean-variable nodes."""
def __init__(self, node_specs=None):
"""Nodes must be ordered with parents before children."""
self.nodes = []
self.variables = []
node_specs = node_specs or []
for node_spec in node_specs:
self.add(node_spec)
def add(self, node_spec):
"""Add a node to the net. Its parents must already be in the
net, and its variable must not."""
node = BayesNode(*node_spec)
assert node.variable not in self.variables
assert all((parent in self.variables) for parent in node.parents)
self.nodes.append(node)
self.variables.append(node.variable)
for parent in node.parents:
self.variable_node(parent).children.append(node)
def variable_node(self, var):
"""Return the node for the variable named var.
>>> burglary.variable_node('Burglary').variable
'Burglary'"""
for n in self.nodes:
if n.variable == var:
return n
raise Exception("No such variable: {}".format(var))
def variable_values(self, var):
"""Return the domain of var."""
return [True, False]
def __repr__(self):
return 'BayesNet({0!r})'.format(self.nodes)
class DecisionNetwork(BayesNet):
"""An abstract class for a decision network as a wrapper for a BayesNet.
Represents an agent's current state, its possible actions, reachable states
and utilities of those states."""
def __init__(self, action, infer):
"""action: a single action node
infer: the preferred method to carry out inference on the given BayesNet"""
super(DecisionNetwork, self).__init__()
self.action = action
self.infer = infer
def best_action(self):
"""Return the best action in the network"""
return self.action
def get_utility(self, action, state):
"""Return the utility for a particular action and state in the network"""
raise NotImplementedError
def get_expected_utility(self, action, evidence):
"""Compute the expected utility given an action and evidence"""
u = 0.0
prob_dist = self.infer(action, evidence, self).prob
for item, _ in prob_dist.items():
u += prob_dist[item] * self.get_utility(action, item)
return u
class InformationGatheringAgent(Agent):
"""
[Figure 16.9]
A simple information gathering agent. The agent works by repeatedly selecting
the observation with the highest information value, until the cost of the next
observation is greater than its expected benefit."""
def __init__(self, decnet, infer, initial_evidence=None):
"""decnet: a decision network
infer: the preferred method to carry out inference on the given decision network
initial_evidence: initial evidence"""
self.decnet = decnet
self.infer = infer
self.observation = initial_evidence or []
self.variables = self.decnet.nodes
def integrate_percept(self, percept):
"""Integrate the given percept into the decision network"""
raise NotImplementedError
def execute(self, percept):
"""Execute the information gathering algorithm"""
self.observation = self.integrate_percept(percept)
vpis = self.vpi_cost_ratio(self.variables)
j = max(vpis)
variable = self.variables[j]
if self.vpi(variable) > self.cost(variable):
return self.request(variable)
return self.decnet.best_action()
def request(self, variable):
"""Return the value of the given random variable as the next percept"""
raise NotImplementedError
def cost(self, var):
"""Return the cost of obtaining evidence through tests, consultants or questions"""
raise NotImplementedError
def vpi_cost_ratio(self, variables):
"""Return the VPI to cost ratio for the given variables"""
v_by_c = []
for var in variables:
v_by_c.append(self.vpi(var) / self.cost(var))
return v_by_c
def vpi(self, variable):
"""Return VPI for a given variable"""
vpi = 0.0
prob_dist = self.infer(variable, self.observation, self.decnet).prob
for item, _ in prob_dist.items():
post_prob = prob_dist[item]
new_observation = list(self.observation)
new_observation.append(item)
expected_utility = self.decnet.get_expected_utility(variable, new_observation)
vpi += post_prob * expected_utility
vpi -= self.decnet.get_expected_utility(variable, self.observation)
return vpi
class BayesNode:
"""A conditional probability distribution for a boolean variable,
P(X | parents). Part of a BayesNet."""
def __init__(self, X, parents, cpt):
"""X is a variable name, and parents a sequence of variable
names or a space-separated string. cpt, the conditional
probability table, takes one of these forms:
* A number, the unconditional probability P(X=true). You can
use this form when there are no parents.
* A dict {v: p, ...}, the conditional probability distribution
P(X=true | parent=v) = p. When there's just one parent.
* A dict {(v1, v2, ...): p, ...}, the distribution P(X=true |
parent1=v1, parent2=v2, ...) = p. Each key must have as many
values as there are parents. You can use this form always;
the first two are just conveniences.
In all cases the probability of X being false is left implicit,
since it follows from P(X=true).
>>> X = BayesNode('X', '', 0.2)
>>> Y = BayesNode('Y', 'P', {T: 0.2, F: 0.7})
>>> Z = BayesNode('Z', 'P Q',
... {(T, T): 0.2, (T, F): 0.3, (F, T): 0.5, (F, F): 0.7})
"""
if isinstance(parents, str):
parents = parents.split()
# We store the table always in the third form above.
if isinstance(cpt, (float, int)): # no parents, 0-tuple
cpt = {(): cpt}
elif isinstance(cpt, dict):
# one parent, 1-tuple
if cpt and isinstance(list(cpt.keys())[0], bool):
cpt = {(v,): p for v, p in cpt.items()}
assert isinstance(cpt, dict)
for vs, p in cpt.items():
assert isinstance(vs, tuple) and len(vs) == len(parents)
assert all(isinstance(v, bool) for v in vs)
assert 0 0:
i_max = backtracking_graph[i - 1][i_max]
return ml_path, ml_probabilities
# _________________________________________________________________________
def fixed_lag_smoothing(e_t, HMM, d, ev, t):
"""
[Figure 15.6]
Smoothing algorithm with a fixed time lag of 'd' steps.
Online algorithm that outputs the new smoothed estimate if observation
for new time step is given.
"""
ev.insert(0, None)
T_model = HMM.transition_model
f = HMM.prior
B = [[1, 0], [0, 1]]
O_t = np.diag(HMM.sensor_dist(e_t))
if t > d:
f = forward(HMM, f, e_t)
O_tmd = np.diag(HMM.sensor_dist(ev[t - d]))
B = matrix_multiplication(np.linalg.inv(O_tmd), np.linalg.inv(T_model), B, T_model, O_t)
else:
B = matrix_multiplication(B, T_model, O_t)
t += 1
if t > d:
# always returns a 1x2 matrix
return [normalize(i) for i in matrix_multiplication([f], B)][0]
else:
return None
# _________________________________________________________________________
def particle_filtering(e, N, HMM):
"""Particle filtering considering two states variables."""
dist = [0.5, 0.5]
# Weight Initialization
w = [0 for _ in range(N)]
# STEP 1
# Propagate one step using transition model given prior state
dist = vector_add(scalar_vector_product(dist[0], HMM.transition_model[0]),
scalar_vector_product(dist[1], HMM.transition_model[1]))
# Assign state according to probability
s = ['A' if probability(dist[0]) else 'B' for _ in range(N)]
w_tot = 0
# Calculate importance weight given evidence e
for i in range(N):
if s[i] == 'A':
# P(U|A)*P(A)
w_i = HMM.sensor_dist(e)[0] * dist[0]
if s[i] == 'B':
# P(U|B)*P(B)
w_i = HMM.sensor_dist(e)[1] * dist[1]
w[i] = w_i
w_tot += w_i
# Normalize all the weights
for i in range(N):
w[i] = w[i] / w_tot
# Limit weights to 4 digits
for i in range(N):
w[i] = float("{0:.4f}".format(w[i]))
# STEP 2
s = weighted_sample_with_replacement(N, s, w)
return s
# _________________________________________________________________________
# TODO: Implement continuous map for MonteCarlo similar to Fig25.10 from the book
class MCLmap:
"""Map which provides probability distributions and sensor readings.
Consists of discrete cells which are either an obstacle or empty"""
def __init__(self, m):
self.m = m
self.nrows = len(m)
self.ncols = len(m[0])
# list of empty spaces in the map
self.empty = [(i, j) for i in range(self.nrows) for j in range(self.ncols) if not m[i][j]]
def sample(self):
"""Returns a random kinematic state possible in the map"""
pos = random.choice(self.empty)
# 0N 1E 2S 3W
orient = random.choice(range(4))
kin_state = pos + (orient,)
return kin_state
def ray_cast(self, sensor_num, kin_state):
"""Returns distance to nearest obstacle or map boundary in the direction of sensor"""
pos = kin_state[:2]
orient = kin_state[2]
# sensor layout when orientation is 0 (towards North)
# 0
# 3R1
# 2
delta = ((sensor_num % 2 == 0) * (sensor_num - 1), (sensor_num % 2 == 1) * (2 - sensor_num))
# sensor direction changes based on orientation
for _ in range(orient):
delta = (delta[1], -delta[0])
range_count = 0
while 0