/*-------------------------------------------------------------------- File : kldiv.c Contents: compute Kullback-Leibler information divergence for simple probability distribution example Author : Christian Borgelt History : 13.10.1999 file created 06.01.2002 log-likelihood computation added 08.01.2002 estimation from data table added --------------------------------------------------------------------*/ #include #include #include /*-------------------------------------------------------------------- Preprocessor Definitions --------------------------------------------------------------------*/ #define LN_2 0.69314718055994530942 /* ln(2) */ #define log2(x) (log(x)/LN_2) #define TPLCNT (sizeof(tab)/sizeof(TUPLE)) /*-------------------------------------------------------------------- Type Definitions --------------------------------------------------------------------*/ typedef struct { /* --- data tuple --- */ int a, b, c; /* attribute values */ int cnt; /* number of occurrences */ } TUPLE; /* (tuple) */ /*-------------------------------------------------------------------- Constants --------------------------------------------------------------------*/ const TUPLE tab[] = { /* --- data tuples */ { 0, 0, 0, 84 }, { 0, 0, 1, 56 }, { 0, 0, 2, 28 }, { 0, 1, 0, 2 }, { 0, 1, 1, 8 }, { 0, 1, 2, 2 }, { 0, 2, 0, 2 }, { 0, 2, 1, 18 }, { 0, 2, 2, 20 }, { 1, 0, 0, 72 }, { 1, 0, 1, 48 }, { 1, 0, 2, 24 }, { 1, 1, 0, 1 }, { 1, 1, 1, 4 }, { 1, 1, 2, 1 }, { 1, 2, 0, 9 }, { 1, 2, 1, 81 }, { 1, 2, 2, 90 }, { 2, 0, 0, 15 }, { 2, 0, 1, 10 }, { 2, 0, 2, 5 }, { 2, 1, 0, 20 }, { 2, 1, 1, 80 }, { 2, 1, 2, 20 }, { 2, 2, 0, 1 }, { 2, 2, 1, 9 }, { 2, 2, 2, 10 }, { 3, 0, 0, 9 }, { 3, 0, 1, 6 }, { 3, 0, 2, 3 }, { 3, 1, 0, 17 }, { 3, 1, 1, 68 }, { 3, 1, 2, 17 }, { 3, 2, 0, 8 }, { 3, 2, 1, 72 }, { 3, 2, 2, 80 } }; /*-------------------------------------------------------------------- Global Variables --------------------------------------------------------------------*/ double p_A[4]; /* marginal distrib. on attribute A */ double p_B[3]; /* marginal distrib. on attribute B */ double p_C[3]; /* marginal distrib. on attribute C */ double p_AB[4][3]; /* marginal distrib. on AxB */ double p_BC[3][3]; /* marginal distrib. on BxC */ double p_AC[4][3]; /* marginal distrib. on AxC */ double p_ABC[4][3][3]; /* full probability distribution */ double p[4][3][3]; /* graphical model distribution */ /*-------------------------------------------------------------------- Functions --------------------------------------------------------------------*/ double kldiv (void) { /* --- Kullback-Leibler info. div. */ int a, b, c; /* attribute values */ double d = 0; /* K-L information divergence */ for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) d += p_ABC[a][b][c] *log2(p_ABC[a][b][c] / p[a][b][c]); return (d > 0) ? d : 0; } /* kldiv() */ /*------------------------------------------------------------------*/ double llhood (void) { /* --- likelihood of dataset */ int i; /* loop variable */ double llh = 0; /* likelihood */ for (i = TPLCNT; --i >= 0; ) llh += tab[i].cnt *log(p[tab[i].a][tab[i].b][tab[i].c]); return llh; } /* llhood() */ /*------------------------------------------------------------------*/ int main (void) { /* --- main function */ int i; /* loop variable */ int a, b, c; /* attribute values */ int sum = 0; /* sum of numbers of occurrences */ /* --- estimate distribution --- */ for (i = TPLCNT; --i >= 0; ) { p_ABC[tab[i].a][tab[i].b][tab[i].c] += tab[i].cnt; sum += tab[i].cnt; } for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p_ABC[a][b][c] /= sum; /* --- compute marginal distributions --- */ for (a = 4; --a >= 0; ) { for (b = 3; --b >= 0; ) { for (c = 3; --c >= 0; ) p_AB[a][b] += p_ABC[a][b][c]; p_A[a] += p_AB[a][b]; } } for (b = 3; --b >= 0; ) { for (c = 3; --c >= 0; ) { for (a = 4; --a >= 0; ) p_BC[b][c] += p_ABC[a][b][c]; p_B[b] += p_BC[b][c]; } } for (c = 3; --c >= 0; ) { for (a = 4; --a >= 0; ) { for (b = 3; --b >= 0; ) p_AC[a][c] += p_ABC[a][b][c]; p_C[c] += p_AC[a][c]; } } /* --- print distributions --- */ printf("AxB|C=c1: AxB|C=c2: AxB|C=c3:\n"); for (b = 3; --b >= 0; ) { for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_ABC[a][b][0]); printf(" "); for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_ABC[a][b][1]); printf(" "); for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_ABC[a][b][2]); printf("\n"); } printf("\nAxB: BxC: AxC:\n"); for (b = 3; --b >= 0; ) { for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_AB[a][b]); printf(" "); for (c = 0; c < 3; c++) printf("%3.0f ", 1000 *p_BC[b][c]); printf(" "); for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_AC[a][b]); printf("\n"); } printf("\nA: "); for (a = 0; a < 4; a++) printf("%3.0f ", 1000 *p_A[a]); printf("\nB: "); for (b = 0; b < 3; b++) printf("%3.0f ", 1000 *p_B[b]); printf("\nC: "); for (c = 0; c < 3; c++) printf("%3.0f ", 1000 *p_C[c]); printf("\n\n"); /* --- Kullback-Leibler information divergence --- */ printf("model KLdiv likelihood\n"); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_A[a] *p_B[b] *p_C[c]; printf("A B C: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_AB[a][b] *p_C[c]; printf("A-B C: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_A[a] *p_BC[b][c]; printf("A B-C: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_AC[a][c] *p_B[b]; printf("A-C B: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_AB[a][b] *p_BC[b][c] /p_B[b]; printf("A-B-C: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_AB[a][b] *p_AC[a][c] /p_A[a]; printf("B-A-C: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_AC[a][c] *p_BC[b][c] /p_C[c]; printf("A-C-B: %.3f %-4.0f\n", kldiv(), llhood()); for (a = 4; --a >= 0; ) for (b = 3; --b >= 0; ) for (c = 3; --c >= 0; ) p[a][b][c] = p_ABC[a][b][c]; printf("A-B-C-A: %.3f %-4.0f\n", kldiv(), llhood()); return 0; } /* main() */