A whole lot of stuff in here, but it is the main functioning script other that the Distributions one.
Try this link for the toolset http://carleton-u.maps.arcgis.com/home/item.html?id=076912624eaa4043bb46ef5068e12626
'''
crosstab.py
Author
Dan Patterson
Dept of Geography and Environmental Studies
Carleton University, Ottawa, Canada
Dan_Patterson@carleton.ca
Date Jan 2010
Modified June 2013
Purpose:
Tools for creating tabulation and crosstabulations of data.
Kappa and Chi-square statistics can also be calculated.
'''
import math
import Distributions as d
import copy
#--------------------------------------------------------------------
class Tabulate(object):
'''
Tabulate a list of data returning the unique values and a frequency
count in the form of a dictionary
'''
def __init__(self, data):
'''
Requires a list of values
'''
assert isinstance(data, (list, tuple)), "Requires a list or tuple"
xs = []
a_dict = {}
for a_key in data:
if a_key not in xs: xs.append(a_key)
if a_dict.has_key(a_key):
a_dict[a_key] += 1
else:
a_dict[a_key] = 1
xs.sort()
self.classes = sorted(set(xs))
self.N = len(self.classes)
self.freq_dict = a_dict
class Crosstab(object):
'''
A simple class that allows one to read a list of x-y data pairs.
Returns the contingency table for use in the Matrix, Kappa or Chi
classes.
'''
def __init__(self, data):
'''
Requires a list of [x, y] pairs which can represent any two variables.
Creates a dictionary of the combinations found.
Returns the unique x values, y values and the dictionary of containing
the pairings and a count for each combination.
'''
assert isinstance(data, (list, tuple)), "Requires a list of lists"
xs = []; ys = []
pairs_dict = {}
for a_pair in data:
x, y = a_pair
a_key = tuple(a_pair)
if x not in xs: xs.append(x)
if y not in ys: ys.append(y)
if pairs_dict.has_key(a_key):
pairs_dict[a_key] += 1
else:
pairs_dict[a_key] = 1
xs.sort(); ys.sort()
all_classes = sorted(set(xs + ys))
combos = [ [x,y] for x in all_classes for y in all_classes ]
self.xs = xs
self.ys = ys
self.pairs_dict = pairs_dict
self.classes = all_classes
self.N = len(all_classes)
self.combos = combos
def count_combos(self):
'''
A pairs-count dictionary of all possible combinations. Returns
the xs, ys as lists and all possible pairs with their counts as a
dictionary.
'''
xyn_list = [] # x, y count list
all_combos = self.combos
pairs_dict = self.pairs_dict
for a_key in all_combos:
a_key = tuple(a_key)
if pairs_dict.has_key(a_key):
xyn_list.append( [a_key[0], a_key[1], pairs_dict[a_key]] )
else:
xyn_list.append( [a_key[0],a_key[1], 0] )
xyn_list.sort()
return xyn_list
class Matrix(object):
'''
A simple matrix class for statistics involving crosstabulation
m = matrix(rows, cols, value) useage
m = matrix() 3x3 matrix with zeros
m = matrix(2,3) 2x3 matrix with zeros
m = matrix(2,3,5) 2x3 matrix with fives
m = matrix([[1,2],[3,4]] matrix of values, r/c from inputs
m.rows matrix rows
m.cols matrix columns
m.data matrix data
'''
def __init__(self, r=None, c=None, value=None):
'''Create a matrix, given an initial value, otherwise 0.'''
if value == None:
value = 0
if r == None:
self.rows = 3
self.cols = 3
self.data = [ [value for x in range(self.cols)] for y in range(self.rows) ]
elif isinstance(r, int) and c==None:
self.rows = r
self.cols = 1
self.data = [ [value for x in range(self.cols)] for y in range(self.rows) ]
elif isinstance(r, int) and isinstance(c, int):
self.rows = r
self.cols = c
self.data = [ [value for x in range(self.cols)] for y in range(self.rows) ]
else:
self.rows = len(r)
self.cols = len(r[0])
self.data = r
def __add__(self,other):
'''
Add two matrices together or a scalar value to a matrix
'''
r = self.rows
c = self.cols
if isinstance(other, (int, float)):
other = Matrix(r,c,other)
r1 = other.rows
c1 = other.cols
else:
r1 = other.rows
c1 = other.cols
assert (r == r1) and (c == c1), "Matrices are not equal in form"
m = Matrix(r,c,0)
for i in range(r):
for j in range(c):
m.data[i][j] = self.data[i][j] + other.data[i][j]
return m
def __div__(self,other):
'''
scalar division
'''
assert isinstance(other, (int, float)) and other > 0.0, "only scalar division is supported"
r = self.rows
c = self.cols
m = Matrix(r,c,0)
for i in range(r):
for j in range(c):
m.data[i][j] = self.data[i][j] / float(other)
return m
def __mul__(self,other):
'''
matrix multiplication using either scalar values or another matrix
returns the resultant matrix and the sum product
'''
r = self.rows
c = self.cols
if isinstance(other, (int, float)):
other = Matrix(r,c,other)
r1 = other.rows
c1 = other.cols
else:
r1 = other.rows
c1 = other.cols
m = Matrix(r,r1,0)
assert (r == c1) , "Rows in matrix A != columns in matrices B"
sumProd = 0.0
for i in range(r):
for j in range(r1):
for k in range(r1):
m.data[i][j] += self.data[i][k] * other.data[k][j]
sumProd += m.data[i][j]
return [m, sumProd]
def __sub__(self,other):
'''
Subtract two matrices or a scalar value from it
'''
r = self.rows
c = self.cols
if isinstance(other, (int, float)):
other = Matrix(r,c,other)
r1 = other.rows
c1 = other.cols
else:
r1 = other.rows
c1 = other.cols
assert (r == r1) and (c == c1), "Matrices are not equal in form"
m = Matrix(r,c,0)
for i in range(r):
for j in range(c):
m.data[i][j] = self.data[i][j] - other.data[i][j]
return m
def __getitem__(self, row, col= None):
'''get a value at row column'''
if col == None:
return self.data[row]
return self.data[row][col]
def __setitem__(self, row, col, value):
'''set row column entry to value'''
self.data[row][col] = value
def __repr__(self):
'''Prints a matrix'''
s = ""
for a_row in self.data:
s += "\n %s" % a_row
return s
def get_col(self, col):
'''Returns the specified column, ordered from 0'''
r = self.rows
c = self.cols
assert col >= 0 and col <= c, "column out of range"
out_col = [ self.data[i][col] for i in range(r) ]
return out_col
def get_row(self, row):
'''Returns the specified row, ordered from 0'''
assert row >= 0 and row <= self.rows, "row out of range"
out_row = self.data[row]
return out_row
def get_trace(self):
'''Returns the trace of a matrix as a list. Matrix must be square'''
assert self.is_square(), "Square matrix required."
trace_list = [ ]
for i in range(self.rows):
trace_list.append(self.data[i][i])
return trace_list
def identity(self):
'''create an identity matrix'''
assert self.is_square(), "Square matrix required."
m = Matrix(self.rows, self.cols, 0)
for i in range(self.rows):
m.data[i][i] = 1
return m
def is_square(self):
'''Test for a square matrix (i.e. 3x3)
useage => assert self.issquare() returns True or False'''
return self.rows == self.cols
def normalize(self, N):
'''Normalizes a square matrix by dividing each entry by its grand total, N
import copy is needed
deepcopy is needed otherwise matrix is altered
norm_matrix = matrix[:] doesn't cut it
'''
norm_matrix = copy.deepcopy(self)
total = float(N)
rows = self.rows; cols = self.cols
for row in range(rows):
for col in range(cols):
norm_matrix.data[row][col] /= total
return norm_matrix
def sum_cols(self):
'''Returns a list of column totals for a matrix.'''
col_totals = []
rows = self.rows; cols = self.cols
for i in range(cols):
temp_list = []
for j in range(rows):
temp_list.append(self.data[j][i])
col_totals.append(sum(temp_list))
return col_totals
def sum_rows(self):
'''Returns a list of row totals for a matrix.'''
row_totals = []
rows = self.rows
for i in range(rows):
row_totals.append(sum(self.data[i]))
return row_totals
def sum_trace(self):
'''
Sum of the diagonal (trace) of a square matrix.
You could simply use: total = sum(get_trace(matrix)) as an alternate'''
assert(self.rows == self.cols), "Non-square matrix"
result = 0
rows = self.rows
for i in range(rows):
result += self.data[i][i]
return result
def sum_matrix(self):
'''Returns the sum of all entries in a matrix'''
total = 0
for row in range(self.rows):
total += sum(self.data[row])
return total
def update_matrix(self, row_col_val_list):
'''
Updates a general matrix with val from the sorted list of
row, col, val numbers as a list of lists of [row,col,val]
'''
n = len(row_col_val_list)
for i in range(n):
r, c, a_val = row_col_val_list[i]
print "r, c, a_val", r, c, a_val
if a_val != 0:
self.data[r - 1][c - 1] = a_val
return self
def update_square(self, row_col_val_list):
'''
Updates a square matrix with val from the sorted list of
row, col, val numbers as a list of lists of [row,col,val]
'''
assert self.is_square(), "Square matrix required"
m = 0
n = self.rows #len(self.data)
for i in range(n):
for j in range(n):
a_val = row_col_val_list[m][2] #[#num = data.count([all_classes[i], all_classes[j]])
self.data[i][j] = a_val
m = m + 1
return self
class Kappa(Matrix):
'''
Matrix functions specific to Kappa statistics.
Inherits from the Matrix class
'''
def kappa_accuracy(self):
'''
Kappa accuracy calculations for a square matrix of counts.
Returns overall, producer and user accuracy
'''
assert self.is_square(), "Square matrix required"
row_totals = self.sum_rows()
col_totals = self.sum_cols()
N = self.sum_matrix()
assert N > 0, "Population size is zero"
trace_list = self.get_trace()
trace_total = sum(trace_list)
overall = trace_total / float(N) #overall accuracy
producer = []
user = []
#overt listing so the logic can be followed
for i in range(len(trace_list)):
if col_totals[i] > 0.0:
producer.append(trace_list[i] / float(col_totals[i]))
else:
producer.append(0.0)
if row_totals[i] > 0.0:
user.append(trace_list[i] / float(row_totals[i]))
else:
user.append(0.0)
return [overall, producer, user]
def kappa_stats(self):
'''
Calculates kappa statics for a normalized square matrix.
Requires the matrix trace and lists of row and column totals.
Returns, kappa, kappa variance and kappa z-score
'''
assert self.is_square(), "Square matrix required"
total = self.sum_matrix()
col_totals = self.sum_cols()
row_totals = self.sum_rows()
trace_list = self.get_trace()
trace_total = self.sum_trace()
N = len(row_totals)
observed = trace_total / float(total) #observed agreement
ii = [ (row_totals[i] * col_totals[i]) for i in range(N) ]
sum_ii = sum( ii )
#kappa = (observed - sum_ii) / (1.0 - sum_ii) #form if normalized
kappa = (total * trace_total - sum_ii) / float(total**2 - sum_ii)
t1 = trace_total / float(total)
t2 = sum_ii / float(total**2)
t3 = [ trace_list[i] * (row_totals[i] + col_totals[i]) for i in range(N) ]
t3 = sum(t3) / float(total**2)
t4 = 0.0
a = 0.0; b = 0.0; c = 0.0; d = 0.0
for i in range(N):
for j in range(N):
if j == i:
a += self.data[i][j]
elif j > i:
b += self.data[i][j]
elif j < i:
c += self.data[i][j]
t4 = t4 + (self.data[i][j] * ((row_totals[i] + col_totals[j])**2))
t4 = t4 / float(total**3)
d = b + c
kappa_var = (t1 * (1.0 - t1)) / ((1.0 - t2)**2) + \
2.0 * (1.0 - t1) * (((2.0 * t1 * t2) - t3)) / ((1.0 - t2)**3) + \
((1.0 - t1)**2) * (t4 - (4.0 * t2**2)) / ((1.0 - t2)**4)
kappa_var = kappa_var / float(total)
kappa_zscore = 0.0
if kappa_var > 0.0:
kappa_zscore = kappa / math.sqrt(kappa_var)
return [kappa, kappa_var, kappa_zscore]
def kappa_confusion(self):
'''
produce a confusion matrix which reflect the possible ways that
observations have been observed and classified
observed
predicted a b a - agreement between observed and predicted
c d b - classed as x, observed not to be x
c - observed as x, not classified as x
d - not classed as x, not observed to be x
sens - sensitivity, equal to producer accuracy
spec - specificity
fpr - falsepositive rate = commission error = 1-specificity
fnr - false negative rate = omission error = 1-sensitivity
ppp - positive predictive power = equal to user accuracy
npp - negative predictive power
conf_matrix - confusion matrix for each class
overall_weighted - sum of conf_matrix for each class
'''
def calc_stats(a,b,c,d):
'''calculate the statistics given the 2x2 matrix of a,b,c,d'''
sens = 0; specif = 0; fpr = 0; fnr = 0; ppp = 0; npp = 0
if (a + c) > 0:
sens = a / float(a + c)
fnr = c / float(a + c)
if (b + d) > 0:
specif = d / float(b + d)
fpr = b / float(b + d)
if (a + b) > 0: ppp = a / float(a + b)
if (c + d) > 0: npp = d / float(c + d)
conf_matrix = Matrix( [ [a, b], [c, d] ] )
return [conf_matrix, sens, specif, fpr, fnr, ppp, npp]
assert self.is_square(), "Square matrix required"
#calculate the stats for each class
N = self.sum_matrix()
trace_list = self.get_trace()
row_totals = self.sum_rows()
col_totals = self.sum_cols()
return_list = []
overall_weighted = Matrix( [[0,0],[0,0]] )
for i in range(len(trace_list)):
a = trace_list[i]
b = row_totals[i] - a
c = col_totals[i] - a
d = N - (a + b + c)
the_returned = calc_stats(a,b,c,d)
conf_matrix, sens, specif, fpr, fnr, ppp, npp = the_returned
overall_weighted = overall_weighted + conf_matrix
return_list.append( [conf_matrix, sens, specif, fpr, fnr, ppp, npp] )
abcd = []
#complete for the overall matrix
for a_row in range(overall_weighted.rows):
for a_col in range(overall_weighted.cols):
abcd.append(overall_weighted.data[a_row][a_col])
a,b,c,d = abcd
the_returned = calc_stats(a,b,c,d)
return_list.append( the_returned )
return return_list
class Chi(Matrix):
'''
Matrix class for chi-square statistics
'''
def chi_expected(self):
'''Calculates the expected chi-squared matrix given the observed matrix'''
rows = self.rows
cols = self.cols
row_totals = self.sum_rows()
col_totals = self.sum_cols()
N = float(sum(row_totals))
exp_matrix = Chi(rows,cols,0.0)
for i in range(rows):
for j in range(cols):
exp_matrix.data[i][j] = (row_totals[i] * col_totals[j]) / N
return exp_matrix
def chi_stats(self, exp_matrix=None):
'''
chi-squared statistic given the observed and expected matrix
returns chi_square value and degrees of freedom
'''
if exp_matrix == None:
exp_matrix = self.chi_expected()
assert self.rows == exp_matrix.rows, "Matrix rows are not equal"
assert self.cols == exp_matrix.cols, "Matrix columns are not equal"
rows = self.rows
cols = self.cols
chi_val = 0.0
O = self.data
E = exp_matrix.data
for i in range(rows):
for j in range(cols):
if E[i][j] != 0:
chi_val += (O[i][j] - E[i][j])**2 / float(E[i][j])
if rows == 1:
df = cols - 1
else:
df = (rows - 1) * (cols - 1)
return [chi_val, df]
#--------------------------------------------------------------------
def kappa_results(data):
'''perform a full kappa result test and return as string'''
header0 = "\n Overall %s \n Producer %s \n User %s"
header1 = "\n kappa %s \n variance %s \n z-score %s"
ct = Crosstab(data)
n = ct.N
xs = ct.xs
ys = ct.ys
xyn_list = ct.count_combos()
sqr_matrix = Kappa(n, n, 0)
matrix = sqr_matrix.update_square(xyn_list)
overall, producer, user = matrix.kappa_accuracy()
kappa, kappa_var, kappa_zscore = matrix.kappa_stats()
results = matrix.kappa_confusion()
class_results = results[0:len(results) - 1] #-1]
overall_weighted = results[len(results) -1 :len(results)]
s = "\n\n............Kappa analysis............\n"
s += "\nClasses information"
s += class_key_results(xs, ys)
s += "\nInput matrix %s\n" % (matrix)
s += "\nAccuracy assessment"
s += header0 % (overall, producer, user)
s += "\n\nKappa statistics\n"
s += header1 % (kappa, kappa_var, kappa_zscore)
s += "\n\nConfusion results\n"
s += confusion_results(class_results)
s += confusion_results(overall_weighted)
s += crosstab_results(matrix)
return s
def chi_results(data):
'''perform a full chi-square test and return the results as a string'''
ct = Crosstab(data)
n = ct.N
xs = ct.xs
ys = ct.ys
xyn_list = ct.count_combos()
sqr_matrix = Chi(n, n, 0)
observed = sqr_matrix.update_square(xyn_list)
expected = observed.chi_expected()
chi_sq, df = observed.chi_stats(expected)
chi_sq = round(chi_sq, 3)
p_value = d.chisq_q(chi_sq, df)
p_value = round(p_value, 4)
s = "\n...........Chi-squared results...........\n"
s += "\nobserved matrix\n" + str(observed)
s += "\n\nexpected matrix\n" + str(expected)
s += "\n\nChi-squared = %s, df = %s, p_val = %s " % (chi_sq, df, p_value)
return s
def confusion_results(class_results):
'''print the results from the confusion matrix'''
headers = ["Confusion matrix:",
"Producer accuracy",
"Specificity",
"Commission error ",
"Omission error ",
"Positive predictive power",
"Negative predictive power"
]
trailers = ["", "Sensitivity", "",
"False Positive Rate = 1 - Specificity",
"False Negative Rate = 1 - Sensitivity",
"User accuracy", ""]
s = ""
j = 0
for a_row in class_results:
if len(class_results) > 1:
s += "\nClass: " + str(j)
else:
s += "\nOverall statistics:"
s += "\n%-25s %10s %-40s" % (headers[0], str(a_row[0]), trailers[0])
for i in range(1,len(a_row)):
a_val = a_row[i]
s += "\n%-25s : %10.6f %-40s" % (headers[i], a_val, trailers[i])
s += "\n"
j += 1
return s
def crosstab_results(matrix, row_totals=None, col_totals=None):
'''
Returns a string allowing you to prints a matrix, a Kappa matrix
or Chi-square matrix or transfer the information to a file.
'''
n = matrix.rows
m = matrix.cols
N = matrix.sum_matrix()
headers = ["Class"]
format_string = "\n%10s "
for i in range(0,n):
headers.append(str(i))
format_string += "%10s "
headers.append("Total")
format_string += "%10s "
if row_totals == None:
row_totals = matrix.sum_rows()
if col_totals == None:
col_totals = matrix.sum_cols()
s = "\nCross-tabulation matrix, observed in columns, classified in rows. \n"
s += format_string % tuple(headers)
for i in range(0,n):
a_row = matrix.data[i]
a_row.insert(0,str(i))
a_row.append(row_totals[i])
s += format_string % tuple(a_row)
last = ["Total"]
last.extend(col_totals)
last.append(N)
s += format_string % tuple(last)
return s
def class_key_results(xs, ys=None):
'''
Prints the class keys to a matrix. Since text may be the data
source, the text are indexed to a number for all further output.
This print statement simply allows one to print the column/row keys.
Requires a list of lists of data pairs
'''
if ys == None:
all_classes = sorted(xs)
else:
all_classes = sorted(set(xs + ys))
s = "\n%10s %10s" % ("Class", "New class")
for i in range(len(all_classes)):
s += "\n%10s %10s" % (all_classes[i], i)
return s
#--------------------------------------------------------------------
if __name__ == "__main__":
'''observed matrix
[0, 3, 1, 3]
[1, 2, 2, 2]
[4, 1, 1, 1]
[1, 1, 1, 2]
'''
data = [[1,2],[1,2],[1,2],[1,3],[1,4],[1,4],[1,4],
[2,1],[2,2],[2,2],[2,3],[2,3],[2,4],[2,4],
[3,1],[3,1],[3,1],[3,1],[3,2],[3,3],[3,4],
[4,1],[4,2],[4,3],[4,4],[4,4]]
s = kappa_results(data)
print s
s = chi_results(data)
print s
No comments:
Post a Comment