diff --git a/test_benchmark.py b/test_benchmark.py index 5287265..bd1270c 100644 --- a/test_benchmark.py +++ b/test_benchmark.py @@ -11,7 +11,8 @@ from cisc import cisc from dc import dc -from test_anm import dc_compat +from dr import dr +from utils import dc_compat, plot_multiline, progress, reverse_argsort from discretizer import * @@ -67,15 +68,19 @@ def load_tubingen_pairs(): def test_tubingen_pairs(): epsilon = 0.0 + level = 0.05 truths = get_ground_truths_of_tubingen_pairs() multivariate_pairs = [52, 53, 54, 55, 71] num_pairs = len(truths) - len(multivariate_pairs) num_correct = 0 num_wrong = 0 + nsample = 0 num_indecisive = 0 + res_cisc, res_dc, res_dr = [], [], [] + diffs_cisc, diffs_dc, diffs_dr = [], [], [] - print "#num\tfound\ttruth\tdelta" + progress(0, 95) for i, data in enumerate(load_tubingen_pairs()): if i + 1 in multivariate_pairs: continue @@ -85,32 +90,91 @@ def test_tubingen_pairs(): # discretizer = UnivariateIPDiscretizer(X, Y) # aX, Xd, aY, Yd = discretizer.discretize() # cisc_score = cisc(Xd, Yd) + nsample += 1 cisc_score = cisc(X, Y) - # cisc_score = dc(dc_compat(Xd), dc_compat(Yd)) - delta = abs(cisc_score[0] - cisc_score[1]) + dc_score = dc(dc_compat(X), dc_compat(Y)) + dr_score = dr(X.tolist(), Y.tolist(), level) + + diffs_cisc.append(abs(cisc_score[0] - cisc_score[1])) + diffs_dc.append(abs(dc_score[0] - dc_score[1])) + diffs_dr.append(abs(dr_score[0] - dr_score[1])) if cisc_score[0] < cisc_score[1]: - cause = "X" - found = "X → Y" + cause_cisc = "X" elif cisc_score[0] > cisc_score[1]: - cause = "Y" - found = "Y → X" + cause_cisc = "Y" + else: + cause_cisc = "" + + if dc_score[0] > dc_score[1]: + cause_dc = "X" + elif dc_score[0] < dc_score[1]: + cause_dc = "Y" + else: + cause_dc = "" + + if dr_score[0] > level and dr_score[1] < level: + cause_dr = "X" + elif dr_score[0] < level and dr_score[1] > level: + cause_dr = "Y" else: - cause = "" - found = "X ~ Y" - - truth_cause = truths[i][0] - truth = "X → Y" if truth_cause == "X" else "Y → X" - if cause == "": - num_indecisive += 1 - elif cause == truth_cause: - num_correct += 1 + cause_dr = "" + + true_cause = truths[i][0] + if cause_cisc == "": + res_cisc.append(random.choice([True, False])) + elif cause_cisc == true_cause: + res_cisc.append(True) else: - num_wrong += 1 + res_cisc.append(False) - print "%4d\t%s\t%s\t%.4f" % (i + 1, found, truth, delta) + if cause_dc == "": + res_dc.append(random.choice([True, False])) + elif cause_dc == true_cause: + res_dc.append(True) + else: + res_dc.append(False) - print "✓ = %3d ✗ = %3d ~ = %3d" % (num_correct, num_wrong, num_indecisive) + if cause_dr == "": + res_dr.append(random.choice([True, False])) + elif cause_dr == true_cause: + res_dr.append(True) + else: + res_dr.append(False) + + progress(nsample, 95) + + # print "✓ = %3d ✗ = %3d ~ = %3d" % (num_correct, num_wrong, + # num_indecisive) + indices_cisc = reverse_argsort(diffs_cisc) + indices_dc = reverse_argsort(diffs_dc) + indices_dr = reverse_argsort(diffs_dr) + + diffs_cisc = [diffs_cisc[i] for i in indices_cisc] + diffs_dc = [diffs_dc[i] for i in indices_dc] + diffs_dr = [diffs_dr[i] for i in indices_dr] + + res_cisc = [res_cisc[i] for i in indices_cisc] + res_dc = [res_dc[i] for i in indices_dc] + res_dr = [res_dr[i] for i in indices_dr] + + dec_rate = np.arange(0.02, 1.01, 0.01) + accs_cisc, accs_dc, accs_dr = [], [], [] + fp = open("results/dec_rate_benchmark.dat", "w") + for r in dec_rate: + maxIdx = int(r * nsample) + rcisc = res_cisc[:maxIdx] + rdc = res_dc[:maxIdx] + rdr = res_dr[:maxIdx] + + accs_cisc.append(sum(rcisc) / len(rcisc)) + accs_dc.append(sum(rdc) / len(rdc)) + accs_dr.append(sum(rdr) / len(rdr)) + fp.write("%.2f %.2f %.2f %.2f" % (r, sum(rdc) / len(rdc), + sum(rdr) / len(rdr), sum(rcisc) / len(rcisc))) + fp.close() + plot_multiline([accs_cisc, accs_dc, accs_dr], dec_rate, [ + "CISC", "DC", "DR"], "decision rate", "accuracy", "decision rate versus accuracy") if __name__ == "__main__": diff --git a/test_synthetic.py b/test_synthetic.py index b758328..9fdc38a 100644 --- a/test_synthetic.py +++ b/test_synthetic.py @@ -4,17 +4,18 @@ """Here, we test the accuracy of various causal inference techniques in identifying ANMs with different ranges of noise term. """ from __future__ import division +from functools import partial import random import sys +from timeit import Timer -import matplotlib as mp -import matplotlib.pyplot as plt import numpy as np from nonparametric_tests import friedman_test, nemenyi_multitest, bonferroni_dunn_test from cisc import cisc from dc import dc from dr import dr +from utils import progress, plot_multiline, reverse_argsort, dc_compat __author__ = "Kailash Budhathoki" @@ -23,16 +24,6 @@ __license__ = "MIT" -def _configure_matplotlib(): - plt.style.use('ggplot') - font_size = 13 - fig = mp.pyplot.gcf() - fig.set_size_inches(18, 12) - - -_configure_matplotlib() - - class InfRes(object): def __init__(self): @@ -45,44 +36,6 @@ def to_str(self, nsample): return str.center("%d %d %d %.2f" % (self.ncorrect, self.nwrong, self.nundec, self.ncorrect / nsample), 20) -def progress(count, total, suffix=''): - bar_len = 60 - filled_len = int(round(bar_len * count / float(total))) - - percents = round(100.0 * count / float(total), 1) - bar = '=' * filled_len + '-' * (bar_len - filled_len) - - sys.stdout.write('[%s] %s%s ...%s\r' % (bar, percents, '%', suffix)) - sys.stdout.flush() # As suggested by Rom Ruben - - -def plot_multiline(Ys, X, labels, xlabel, ylabel, title, fig_name=None): - plt.title(title) - plt.xlabel(xlabel) - plt.ylabel(ylabel) - - for i, Y in enumerate(Ys): - plt.plot(X, Y, label=labels[i], marker="o", linewidth=2.0, - dash_joinstyle="bevel", solid_capstyle="round", markeredgecolor="none", linestyle="-") - plt.legend(loc="upper left") - plt.ylim(0.0, 1.001) - if fig_name: - plt.savefig("results/%s" % fig_name) - else: - plt.show() - plt.cla() - - -def reverse_argsort(X): - indices = range(len(X)) - indices.sort(key=X.__getitem__, reverse=True) - return indices - - -def dc_compat(X): - return [[x] for x in X] - - def map_randomly(Xd, fXd): f = dict() for x in Xd: @@ -91,24 +44,24 @@ def map_randomly(Xd, fXd): return f -def generate_variable(src, size): - if src == "unif": +def generate_X(src, size): + if src == "uniform": maxX = random.randint(1, 10) X = [random.randint(1, maxX) for i in xrange(size)] - elif src == "multinom": + elif src == "multinomial": p_nums = [random.randint(1, 10) for i in xrange(random.randint(1, 11))] p_vals = [v / sum(p_nums) for v in p_nums] X = np.random.multinomial(size, p_vals, size=1)[0].tolist() X = [[i + 1] * f for i, f in enumerate(X)] X = [j for sublist in X for j in sublist] - elif src == "binom": + elif src == "binomial": n = random.randint(1, 40) p = random.uniform(0.1, 0.9) X = np.random.binomial(n, p, size).tolist() - elif src == "geom": + elif src == "geometric": p = random.uniform(0.1, 0.9) X = np.random.geometric(p, size).tolist() - elif src == "hypgeom": + elif src == "hypergeometric": ngood = random.randint(1, 41) nbad = random.randint(1, 41) nsample = random.randint(1, min(41, ngood + nbad)) @@ -116,31 +69,33 @@ def generate_variable(src, size): elif src == "poisson": lam = random.randint(1, 10) X = np.random.poisson(lam, size).tolist() - elif src == "nbinom": + elif src == "negativeBinomial": n = random.randint(1, 40) p = random.uniform(0.1, 0.9) X = np.random.negative_binomial(n, p, size).tolist() return X -# -------------------------- all the tests are below -------------------------- +def generate_additive_N(size): + t = random.randint(1, 7) + suppN = range(-t, t + 1) + N = [random.choice(suppN) for i in xrange(size)] + return N def test_decision_rate(): - srcs = ["unif", "multinom", "binom", - "geom", "hypgeom", "poisson", "nbinom"] + srcsX = ["uniform", "binomial", "negativeBinomial", + "geometric", "hypergeometric", "poisson", "multinomial"] total = len(srcs) * len(srcs) count = 0 progress(count, total) - for srcX in srcs: - for srcN in srcs: - count += 1 - _test_decision_rate(srcX, srcN) - progress(count, total) + for srcX in srcsX: + _decision_rate(srcX) + progress(count, total) -def _test_decision_rate(srcX, srcN): - nsample = 500 +def _decision_rate(srcX): + nsample = 1000 size = 1000 level = 0.05 suppfX = range(-7, 8) @@ -149,10 +104,10 @@ def _test_decision_rate(srcX, srcN): res_cisc, res_dc, res_dr = [], [], [] for k in xrange(nsample): - X = generate_variable(srcX, size) + X = generate_X(srcX, size) suppX = list(set(X)) f = map_randomly(suppX, suppfX) - N = generate_variable(srcN, size) + N = generate_additive_N(size) Y = [f[X[i]] + N[i] for i in xrange(size)] cisc_score = cisc(X, Y) @@ -210,70 +165,79 @@ def _test_decision_rate(srcX, srcN): acc_cisc.append(acc_cisc_r) acc_dc.append(acc_dc_r) acc_dr.append(acc_dr_r) - + print dec_rate + print acc_dc + print acc_dr + print acc_cisc + + with open("results/dec_rate_synth_data.dat", "w") as writer: + for i, r in enumerate(dec_rate): + writer.write("%.2f %.2f %.2f %.2f\n" % + (r, acc_dc[i], acc_dr[i], acc_cisc[i])) plot_multiline([acc_cisc, acc_dc, acc_dr], dec_rate, [ - "CISC", "DC", "DR"], "decision rate", "accuracy", "decision rate versus accuracy", "%sX-%sN.png" % (srcX, srcN)) + "CISC", "DC", "DR"], "decision rate", "accuracy", "decision rate versus accuracy", "dec_rate_%sX.png" % srcX) -def test_anm(): - nsim = 500 +def test_accuracy(): + nsim = 1000 size = 1000 level = 0.05 suppfX = range(-7, 8) - srcs = ["unif", "multinom", "binom", - "geom", "hypgeom", "poisson", "nbinom"] - - for srcX in srcs: - print "typeX :: %s" % srcX - print "%10s%10s%10s%10s" % ("typeN", "acc(cisc)", "acc(dc)", "acc(dr)") - - for srcN in srcs: # for each type of N - ncorrect_cisc, ncorrect_dc, ncorrect_dr = 0, 0, 0 - - for k in xrange(nsim): - X = generate_variable(srcX, size) - suppX = list(set(X)) - f = map_randomly(suppX, suppfX) - N = generate_variable(srcN, size) - Y = [f[X[i]] + N[i] for i in xrange(size)] + srcsX = ["uniform", "binomial", "negativeBinomial", + "geometric", "hypergeometric", "poisson", "multinomial"] + print "-" * 48 + print "%18s%10s%10s%10s" % ("TYPE_X", "DC", "DR", "CISC") + print "-" * 48 + + fp = open("results/accuracy_synthetic.dat", "w") + for srcX in srcsX: + ncorrect_cisc, ncorrect_dc, ncorrect_dr = 0, 0, 0 + for k in xrange(nsim): + X = generate_X(srcX, size) + suppX = list(set(X)) + f = map_randomly(suppX, suppfX) + N = generate_additive_N(size) + Y = [f[X[i]] + N[i] for i in xrange(size)] - cisc_score = cisc(X, Y) - dc_score = dc(dc_compat(X), dc_compat(Y)) - dr_score = dr(X, Y, level) + cisc_score = cisc(X, Y) + dc_score = dc(dc_compat(X), dc_compat(Y)) + dr_score = dr(X, Y, level) - ncorrect_cisc += int(cisc_score[0] < cisc_score[1]) - ncorrect_dc += int(dc_score[0] > dc_score[1]) - ncorrect_dr += int(dr_score[0] > level and dr_score[1] < level) + ncorrect_cisc += int(cisc_score[0] < cisc_score[1]) + ncorrect_dc += int(dc_score[0] > dc_score[1]) + ncorrect_dr += int(dr_score[0] > level and dr_score[1] < level) - print "%10s%10.2f%10.2f%10.2f" % (srcN, ncorrect_cisc / nsim, ncorrect_dc / nsim, ncorrect_dr / nsim) - print + print "%18s%10.2f%10.2f%10.2f" % (srcX, ncorrect_dc / nsim, ncorrect_dr / nsim, ncorrect_cisc / nsim) + fp.write("%s %.2f %.2f %.2f\n" % (srcX, ncorrect_dc / + nsim, ncorrect_dr / nsim, ncorrect_cisc / nsim)) + print "-" * 48 + fp.close() -def test_sc_vs_anm_space(): +def test_accuracy_at_model_space(): nsim = 500 size = 1500 level = 0.05 - fXd = range(-7, 8) + suppfX = range(-7, 8) S = range(1, 7) - srcs = ["unif", "multinom", "binom", - "geom", "hypgeom", "poisson", "nbinom"] + srcsX = ["uniform", "binomial", "negativeBinomial", + "geometric", "hypergeometric", "poisson", "multinomial"] - print "%s | %s | %s | %s" % ("typeX".center(8), "sc_space".center(62), "dr_space".center(62), "full_space".center(62)) - print "%8s | %s %s %s | %s %s %s | %s %s %s" % ("", "sc".center(20), "dc".center(20), "dr".center(20), "sc".center(20), "dc".center(20), "dr".center(20), "sc".center(20), "dc".center(20), "dr".center(20)) + print "%s | %s | %s | %s" % ("typeX".center(18), "sc_space".center(62), "dr_space".center(62), "full_space".center(62)) + print "%18s | %s %s %s | %s %s %s | %s %s %s" % ("", "sc".center(20), "dc".center(20), "dr".center(20), "sc".center(20), "dc".center(20), "dr".center(20), "sc".center(20), "dc".center(20), "dr".center(20)) - for src in srcs: + for srcX in srcsX: nsample_sc, nsample_anm, nsample_all = 0, 0, 0 rsc_sc, rsc_dc, rsc_dr = InfRes(), InfRes(), InfRes() ranm_sc, ranm_dc, ranm_dr = InfRes(), InfRes(), InfRes() rall_sc, rall_dc, rall_dr = InfRes(), InfRes(), InfRes() for k in xrange(nsim): - X = generate_variable(src, size) - Xd = list(set(X)) - f = map_randomly(Xd, fXd) - t = random.choice(S) - N = range(-t, t + 1) - Y = [f[X[i]] + random.choice(N) for i in xrange(size)] + X = generate_X(srcX, size) + suppX = list(set(X)) + f = map_randomly(suppX, suppfX) + N = generate_additive_N(size) + Y = [f[X[i]] + N[i] for i in xrange(size)] cisc_score = cisc(X, Y) dc_score = dc(dc_compat(X), dc_compat(Y)) @@ -328,19 +292,19 @@ def test_sc_vs_anm_space(): rall_dr.nwrong += int(dr_score[0] < level and dr_score[1] > level) rall_dr.nundec += int(undecided_dr) - print "%s | %s %s %s | %s %s %s | %s %s %s" % (src.center(8), rsc_sc.to_str(nsample_sc), rsc_dc.to_str(nsample_sc), rsc_dr.to_str(nsample_sc), ranm_sc.to_str(nsample_anm), ranm_dc.to_str(nsample_anm), ranm_dr.to_str(nsample_anm), rall_sc.to_str(nsample_all), rall_dc.to_str(nsample_all), rall_dr.to_str(nsample_all)) + print "%s | %s %s %s | %s %s %s | %s %s %s" % (srcX.center(18), rsc_sc.to_str(nsample_sc), rsc_dc.to_str(nsample_sc), rsc_dr.to_str(nsample_sc), ranm_sc.to_str(nsample_anm), ranm_dc.to_str(nsample_anm), ranm_dr.to_str(nsample_anm), rall_sc.to_str(nsample_all), rall_dc.to_str(nsample_all), rall_dr.to_str(nsample_all)) -def test_anm_nonparam_tests(): +def test_accuracy_with_nonparam_tests(): nsim = 250 size = 1000 level = 0.05 fXd = range(-7, 8) S = range(1, 7) - srcs = ["unif", "multinom", "binom", - "geom", "hypgeom", "poisson", "nbinom"] + srcs = ["uniform", "binomial", "negativeBinomial", + "geometric", "hypergeometric", "poisson", "multinomial"] - print "%10s%10s%10s%10s" % ("X_TYPE", "ACC(CISC)", "ACC(DC)", "ACC(DR)") + print "%18s%10s%10s%10s" % ("X_TYPE", "ACC(CISC)", "ACC(DC)", "ACC(DR)") for src in srcs: ncorrect_this, nwrong_this, nind_this = 0, 0, 0 @@ -349,7 +313,7 @@ def test_anm_nonparam_tests(): diff_cisc, diff_dc, diff_dr = [], [], [] for k in xrange(nsim): - X = generate_variable(src, size) + X = generate_X(src, size) Xd = list(set(X)) f = map_randomly(Xd, fXd) t = random.choice(S) @@ -395,8 +359,90 @@ def test_anm_nonparam_tests(): else: print "all the algorithms are the same" - print "%10s%10.2f%10.2f%10.2f" % (src, ncorrect_this / nsim, ncorrect_dc / nsim, ncorrect_dr / nsim) + print "%18s%10.2f%10.2f%10.2f" % (src, ncorrect_this / nsim, ncorrect_dc / nsim, ncorrect_dr / nsim) + + +def test_runtime(): + nrun = 5 + level = 0.05 + suppfX = range(-7, 8) + nvalues = 10 + sizes = [100000, 200000, 300000, 400000, 500000, + 600000, 700000, 800000, 900000, 1000000] + + for size in sizes: + p_nums = [random.randint(1, 10) for i in xrange(nvalues)] + p_vals = [v / sum(p_nums) for v in p_nums] + X = np.random.multinomial(size, p_vals, size=1)[0].tolist() + X = [[i + 1] * f for i, f in enumerate(X)] + X = [j for sublist in X for j in sublist] + suppX = list(set(X)) + f = map_randomly(suppX, suppfX) + N = generate_additive_N(size) + Y = [f[X[i]] + N[i] for i in xrange(size)] + dcX, dcY = dc_compat(X), dc_compat(Y) + + dc_t = Timer(partial(dc, dcX, dcY), + setup="""from __main__ import *""").timeit(nrun) + dr_t = Timer(partial(dr, X, Y, level), + setup="""from __main__ import *""").timeit(nrun) + cisc_t = Timer(partial(cisc, X, Y), + setup="""from __main__ import *""").timeit(nrun) + print size, dc_t, dr_t, cisc_t + + +def test_power(): + sizes = [1000, 2000, 3000, 4000, 5000] + level = 0.05 + srcX = "multinomial" + nsim_cutoff, nsim_power = 100, 100 + suppfX = range(-7, 8) + + for size in sizes: + diffs_dc, diffs_dr, diffs_cisc = [], [], [] + for i in xrange(nsim_cutoff): + X = generate_X(srcX, size) + suppX = list(set(X)) + f = map_randomly(suppX, suppfX) + Y = [f[X[i]] for i in xrange(size)] + + dc_score = dc(dc_compat(X), dc_compat(Y)) + dr_score = dr(X, Y, level) + cisc_score = cisc(X, Y) + + diffs_dc.append(abs(dc_score[0] - dc_score[1])) + diffs_dr.append(abs(dr_score[0] - dr_score[1])) + diffs_cisc.append(abs(cisc_score[0] - cisc_score[1])) + + cutoff_dc = np.percentile(diffs_dc, 5) + cutoff_dr = np.percentile(diffs_dr, 5) + cutoff_cisc = np.percentile(diffs_cisc, 5) + print cutoff_dc, cutoff_dr, cutoff_cisc + + diffs_dc, diffs_dr, diffs_cisc = [], [], [] + for i in xrange(nsim_power): + X = generate_X(srcX, size) + suppX = list(set(X)) + f = map_randomly(suppX, suppfX) + N = generate_additive_N(size) + Y = [f[X[i]] + N[i] for i in xrange(size)] + + dc_score = dc(dc_compat(X), dc_compat(Y)) + dr_score = dr(X, Y, level) + cisc_score = cisc(X, Y) + + diffs_dc.append(abs(dc_score[0] - dc_score[1])) + diffs_dr.append(abs(dr_score[0] - dr_score[1])) + diffs_cisc.append(abs(cisc_score[0] - cisc_score[1])) + + power_dc = sum(diffs_dc > cutoff_dc) * 1.0 / nsim_power + power_dr = sum(diffs_dr > cutoff_dr) * 1.0 / nsim_power + power_cisc = sum(diffs_cisc > cutoff_cisc) * 1.0 / nsim_power + + print size, power_dc, power_dr, power_cisc if __name__ == "__main__": - test_sc_vs_anm_space() + test_power() + # test_runtime() + # _decision_rate("multinomial") diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..dd9d825 --- /dev/null +++ b/utils.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +"""All the utility functions used everywhere in this project are here. +""" +import sys + +import matplotlib as mp +import matplotlib.pyplot as plt + + +__author__ = "Kailash Budhathoki" +__email__ = "kbudhath@mpi-inf.mpg.de" +__copyright__ = "Copyright (c) 2017" +__license__ = "MIT" + + +def _configure_matplotlib(): + plt.style.use('ggplot') + font_size = 13 + fig = mp.pyplot.gcf() + fig.set_size_inches(18, 12) + + +_configure_matplotlib() + + +def progress(count, total, suffix=''): + bar_len = 60 + filled_len = int(round(bar_len * count / float(total))) + + percents = round(100.0 * count / float(total), 1) + bar = '=' * filled_len + '-' * (bar_len - filled_len) + + sys.stdout.write('[%s] %s%s ...%s\r' % (bar, percents, '%', suffix)) + sys.stdout.flush() # As suggested by Rom Ruben + + +def plot_multiline(Ys, X, labels, xlabel, ylabel, title, fig_name=None): + plt.title(title) + plt.xlabel(xlabel) + plt.ylabel(ylabel) + + for i, Y in enumerate(Ys): + plt.plot(X, Y, label=labels[i], marker="o", linewidth=2.0, + dash_joinstyle="bevel", solid_capstyle="round", markeredgecolor="none", linestyle="-") + plt.legend(loc="upper left") + plt.ylim(0.0, 1.001) + if fig_name: + plt.savefig("results/%s" % fig_name) + else: + plt.show() + plt.cla() + + +def reverse_argsort(X): + indices = range(len(X)) + indices.sort(key=X.__getitem__, reverse=True) + return indices + + +def dc_compat(X): + return [[x] for x in X]