Part of iPyMacLern project.
Copyright (C) 2016 by Eka A. Kurniawan
eka.a.kurniawan(ta)gmail(tod)com
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.
In [1]:
# Display graph inline
%matplotlib inline
# Display graph in 'retina' format for Mac with retina display. Others, use PNG or SVG format.
%config InlineBackend.figure_format = 'retina'
#%config InlineBackend.figure_format = 'PNG'
#%config InlineBackend.figure_format = 'SVG'
In [2]:
import sys
print("Python %s" % sys.version)
In [3]:
import numpy as np
print("NumPy %s" % np.__version__)
In [4]:
import sklearn
from sklearn.utils.extmath import cartesian
print("scikit-learn %s" % sklearn.__version__)
In [5]:
import pandas as pd
print("Pandas %s" % pd.__version__)
In [6]:
import math
Factor $\phi$ is a function mapping values from a set of random variables $D$, $Val(D)$ to real number $I\!R$.$^{[1,2]}$
$$\phi(X_1, \ldots, X_k):Val(X_1, \ldots, X_k) \rightarrow I\!R$$Where $D$ is also called the scope of the factor.
$$D = Scope[\phi] = \{X_1, \ldots, X_k\}$$Data structure for factor $\phi$ with three binary variables $X_2$, $X_1$, and $X_3$ that are all assigned to 1, are written as follow.$^{[3]}$
In [7]:
phi = {'variables': np.array(['X_2', 'X_1', 'X_3']),
'cardinalities': np.array([2, 2, 2]),
'values': np.ones(8)}
The card key stores cardinality of all variables. In the example, the cardinality of variable $X_3$, $|Val(X_3)|$ is 2.
Following function converts factor data structure into the table form for displaying purpose.
In [8]:
def convert_factor_ds_to_table(phi):
cart = cartesian([range(d) for d in phi['cardinalities']])
# Construct table
df = pd.DataFrame(cart)
df.columns = ['$%s$' % v for v in phi['variables']]
df['$\phi(%s)$' % ','.join([v for v in phi['variables']])] = phi['values']
return df
In [9]:
phi_table = convert_factor_ds_to_table(phi)
In [10]:
phi_table
Out[10]:
Convert from factor assignment $A$ to the value index $I$ based on the cardinalities $D$.
In [11]:
def convert_assignment_to_index(A, D):
step = [np.prod(D[i+1:]) for i in range(len(D)-1)] + [0]
return np.sum(np.multiply(A, step)) + A[-1]
In [12]:
A = [1, 0, 1]
I = convert_assignment_to_index(A, D = phi['cardinalities'])
I
Out[12]:
Convert from factor value index $I$ to the assignment $A$ based on the cardinalities $D$.
In [13]:
def convert_index_to_assignment(I, D):
step = [np.prod(D[i+1:]) for i in range(len(D)-1)] + [0]
I_tmp = I
A = []
for i in range(len(D)-1):
a = int(I_tmp / step[i])
A.append(a)
I_tmp = I_tmp - (step[i] * a)
A.append(I_tmp)
return A
In [14]:
I = 5
A = convert_index_to_assignment(I, D = phi['cardinalities'])
A
Out[14]:
In [15]:
class factor:
def __init__(self, variables, cardinalities, values):
self.variables = np.array(variables, '<U255')
self.cardinalities = np.array(cardinalities)
self.values = np.array(values)
def get_table(self):
cart = cartesian([range(cardinality) for cardinality in self.cardinalities])
# Construct table
df = pd.DataFrame(cart)
df.columns = ['$%s$' % l for l in self.variables]
df['$\phi(%s)$' % ','.join([l for l in self.variables])] = self.values
return df
def get_index_from_assignment(self, A):
step = [np.prod(self.cardinalities[i+1:]) for i in range(len(self.cardinalities)-1)] + [0]
return np.sum(np.multiply(A, step)) + A[-1]
def get_assignment_from_index(self, I):
step = [np.prod(self.cardinalities[i+1:]) for i in range(len(self.cardinalities)-1)] + [0]
I_tmp = I
A = []
for i in range(len(self.cardinalities)-1):
a = int(I_tmp / step[i])
A.append(a)
I_tmp = I_tmp - (step[i] * a)
A.append(I_tmp)
return A
def get_value_of_assignment(self, A):
I = self.get_index_from_assignment(A)
return self.values[I]
def set_value_of_assignment(self, A, v):
I = self.get_index_from_assignment(A)
self.values[I] = v
def reduce(self, E):
new_values = self.values.copy()
for variable, value in E.items():
if variable not in self.variables:
continue
variable_index = np.where(self.variables == variable)[0][0]
cart = cartesian([range(cardinality) for cardinality in self.cardinalities])
reduce_index = np.where(cart[:,variable_index] != value)[0]
new_values[reduce_index] = 0
return factor(variables = self.variables,
cardinalities = self.cardinalities,
values = new_values)
In [16]:
phi = factor(variables = ['X_2', 'X_1', 'X_3'],
cardinalities = [2, 2, 2],
values = np.ones(8))
In [17]:
phi.variables
Out[17]:
In [18]:
phi.cardinalities
Out[18]:
In [19]:
phi.values
Out[19]:
In [20]:
phi.get_index_from_assignment([1, 0, 1])
Out[20]:
In [21]:
phi.get_assignment_from_index(5)
Out[21]:
In [22]:
phi.get_value_of_assignment([1, 0, 1])
Out[22]:
In [23]:
phi.set_value_of_assignment([1, 0, 1], 6)
phi.get_value_of_assignment([1, 0, 1])
Out[23]:
get_member_index function returns the lowest indices of values in array $A$ found in array $B$.
In [24]:
def get_member_index(A, B):
C = []
for v in A:
idx = np.where(B == v)[0]
if len(idx):
C.append(idx[0])
return np.array(C)
In [25]:
A = np.array([1,2,3])
B = np.array([2,1,3])
get_member_index(A, B)
Out[25]:
In [26]:
A = np.array([1,2,3])
B = np.array([1,4,3,5,2])
get_member_index(A, B)
Out[26]:
Sample factors ($\phi(X_1)$, $\phi(X_2)$, and $\phi(X_3)$).$^{[3]}$
In [27]:
phi_X1 = factor(variables = ['X_1'],
cardinalities = [2],
values = [0.11, 0.89])
phi_X2 = factor(variables = ['X_1', 'X_2'],
cardinalities = [2, 2],
values = [0.59, 0.41, 0.22, 0.78])
phi_X3 = factor(variables = ['X_2', 'X_3'],
cardinalities = [2, 2],
values = [0.39, 0.61, 0.06, 0.94])
In [28]:
phi_X1.get_table()
Out[28]:
In [29]:
phi_X2.get_table()
Out[29]:
In [30]:
phi_X3.get_table()
Out[30]:
factor_product function returns product of two factors $A$ and $B$ stored in $C$.
In [31]:
def factor_product(A, B):
variables = np.union1d(A.variables, B.variables)
ttl_variables = len(variables)
mapA = get_member_index(A.variables, variables)
mapB = get_member_index(B.variables, variables)
cardinalities = np.zeros(ttl_variables, int)
cardinalities[mapA] = A.cardinalities
cardinalities[mapB] = B.cardinalities
ttl_values = np.prod(cardinalities)
values = np.zeros(ttl_values)
C = factor(variables = variables,
cardinalities = cardinalities,
values = values)
assignments = np.array([C.get_assignment_from_index(i) for i in range(ttl_values)])
indexA = [A.get_index_from_assignment(assignment) for assignment in assignments[:, mapA]];
indexB = [B.get_index_from_assignment(assignment) for assignment in assignments[:, mapB]];
values = np.multiply(A.values[indexA], B.values[indexB])
for i in range(len(C.values)):
C.set_value_of_assignment(assignments[i], values[i])
return C
In [32]:
phi_X1_X2 = factor_product(phi_X1, phi_X2)
In [33]:
phi_X1_X2.get_table()
Out[33]:
factor_marginalization function returns new factor $B$ in which variables in factor $A$ have been summed up based on variables in $V$.
In [34]:
def factor_marginalization(A, V):
variables = np.setdiff1d(A.variables, V)
ttl_variables = len(variables)
mapB = get_member_index(variables, A.variables)
cardinalities = np.zeros(ttl_variables, int)
cardinalities = A.cardinalities[mapB]
ttl_values = np.prod(cardinalities)
values = np.zeros(ttl_values)
B = factor(variables = variables,
cardinalities = cardinalities,
values = values)
assignments = np.array([A.get_assignment_from_index(i) for i in range(len(A.values))])
indics = np.array([B.get_index_from_assignment(assignment[mapB]) for assignment in assignments])
for i in range(len(B.values)):
value = np.sum(A.values[np.where(indics == i)])
B.values[i] = value
return B
In [35]:
phi_X2_marginalized = factor_marginalization(phi_X2, ['X_2'])
In [36]:
phi_X2_marginalized.get_table()
Out[36]:
In [37]:
evidence = {
'X_2': 0,
'X_3': 1,
}
In [38]:
phi_X1.reduce(evidence).get_table()
Out[38]:
In [39]:
phi_X2.reduce(evidence).get_table()
Out[39]:
In [40]:
phi_X3.reduce(evidence).get_table()
Out[40]:
In [41]:
def calculate_joint_distribution(F):
joint = F[0]
for i in range(1, len(F)):
joint = factor_product(joint, F[i])
return joint
In [42]:
F = [phi_X1, phi_X2, phi_X3]
phi_X1_X2_X3 = calculate_joint_distribution(F)
In [43]:
phi_X1_X2_X3.get_table()
Out[43]:
In [44]:
def compare_factors(A, B):
ttl_variables_A = len(A.variables)
map_variables = get_member_index(A.variables, B.variables)
# Compare variables
if np.array_equal(A.variables, B.variables[map_variables]):
print("%-58s%12s" % ("Variables:", "same"))
else:
print("%-58s%12s" % ("Variables:", "DIFFERENT"))
print(" A: ", A.variables)
print(" B: ", B.variables[map_variables])
# Compare cardinalities
if np.array_equal(A.cardinalities, B.cardinalities[map_variables]):
print("%-58s%12s" % ("Cardinalities:", "same"))
else:
print("%-58s%12s" % ("Cardinalities:", "DIFFERENT"))
print(" A: ", A.cardinalities)
print(" B: ", B.cardinalities[map_variables])
# Compare values
ttl_values = np.prod(A.cardinalities)
cart = cartesian([range(cardinality) for cardinality in A.cardinalities])
for i in range(ttl_values):
variables_assignments = ["%s=%s" % (A.variables[j], cart[i][j]) for j in range(ttl_variables_A)]
if math.isclose(A.get_value_of_assignment(cart[i]),
B.get_value_of_assignment(cart[i][map_variables]),
rel_tol=1e-15, abs_tol=0.0):
print("%-58s%12s" % ("Values(" + ",".join(variables_assignments) + "):", "same"))
else:
print("%-58s%12s" % ("Values(" + ",".join(variables_assignments) + "):", "DIFFERENT"))
print(" A: ", A.get_value_of_assignment(cart[i]))
print(" B: ", B.get_value_of_assignment(cart[i][map_variables]))
In [45]:
phi_A_B = factor(variables = ['A', 'B'],
cardinalities = [3, 2],
values = [0.5, 0.8, 0.1, 0.0, 0.3, 0.9])
phi_B_C = factor(variables = ['B', 'C'],
cardinalities = [2, 2],
values = [0.5, 0.7, 0.1, 0.2])
phi_A_B_C_expected1 = factor(variables = ['A', 'B', 'C'],
cardinalities = [3, 2, 2],
values = [0.25, 0.35, 0.08, 0.16, 0.05, 0.07, 0.00, 0.00, 0.15, 0.21, 0.09, 0.18])
phi_A_B_C_expected2 = factor(variables = ['B', 'A', 'C'],
cardinalities = [2, 3, 2],
values = [0.25, 0.35, 0.05, 0.07, 0.15, 0.21, 0.08, 0.16, 0.00, 0.00, 0.09, 0.18])
In [46]:
phi_A_B_C = factor_product(phi_A_B, phi_B_C)
In [47]:
compare_factors(phi_A_B_C, phi_A_B_C_expected1)
In [48]:
compare_factors(phi_A_B_C, phi_A_B_C_expected2)
In [49]:
phi_A_B_C = factor(variables = ['A', 'B', 'C'],
cardinalities = [3, 2, 2],
values = [0.25, 0.35, 0.08, 0.16, 0.05, 0.07, 0.00, 0.00, 0.15, 0.21, 0.09, 0.18])
phi_A_C_expected = factor(variables = ['A', 'C'],
cardinalities = [3, 2],
values = [0.33, 0.51, 0.05, 0.07, 0.24, 0.39])
In [50]:
phi_A_C = factor_marginalization(phi_A_B_C, ['B'])
In [51]:
compare_factors(phi_A_C, phi_A_C_expected)
In [52]:
evidence = {
'C': 0,
}
phi_A_B_C = factor(variables = ['A', 'B', 'C'],
cardinalities = [3, 2, 2],
values = [0.25, 0.35, 0.08, 0.16, 0.05, 0.07, 0.00, 0.00, 0.15, 0.21, 0.09, 0.18])
phi_A_B_C0_expected = factor(variables = ['A', 'B', 'C'],
cardinalities = [3, 2, 2],
values = [0.25, 0.0, 0.08, 0.0, 0.05, 0.0, 0.00, 0.0, 0.15, 0.0, 0.09, 0.0])
In [53]:
phi_A_B_C0 = phi_A_B_C.reduce(evidence)
In [54]:
compare_factors(phi_A_B_C0, phi_A_B_C0_expected)
In [55]:
phi_A_B = factor(variables = ['A', 'B'],
cardinalities = [3, 2],
values = [0.5, 0.8, 0.1, 0.0, 0.3, 0.9])
phi_B_C = factor(variables = ['B', 'C'],
cardinalities = [2, 2],
values = [0.5, 0.7, 0.1, 0.2])
phi_A_B_C_expected = factor(variables = ['A', 'B', 'C'],
cardinalities = [3, 2, 2],
values = [0.25, 0.35, 0.08, 0.16, 0.05, 0.07, 0.00, 0.00, 0.15, 0.21, 0.09, 0.18])
In [56]:
F = [phi_A_B]
phi_A_B_joint = calculate_joint_distribution(F)
In [57]:
compare_factors(phi_A_B_joint, phi_A_B)
In [58]:
F = [phi_A_B, phi_B_C]
phi_A_B_C = calculate_joint_distribution(F)
In [59]:
compare_factors(phi_A_B_C, phi_A_B_C_expected)