condmutualinfo.java

bayes network classifier toolbox 贝叶斯网络分类工具箱

/**
 *  JBNC - Bayesian Network Classifiers Toolbox <p>
 *
 *  Latest release available at http://sourceforge.net/projects/jbnc/ <p>
 *
 *  Copyright (C) 1999-2003 Jarek Sacha <p>
 *
 *  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 2 of the License, or (at your option)
 *  any later version. <p>
 *
 *  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. <p>
 *
 *  You should have received a copy of the GNU General Public License along with
 *  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
 *  Place - Suite 330, Boston, MA 02111-1307, USA. <br>
 *  http://www.fsf.org/licenses/gpl.txt
 */
package jbnc.util;

import jbnc.dataset.AttributeSpecs;

/**
 *  Estimate a joint probability distribution of three discreate variables
 *
 * @author     Jarek Sacha
 * @since      June 1, 1999
 */
public class CondMutualInfo {
    private final static double inv_log2 = 1.0 / Math.log(2);

    private AttributeSpecs[] names;
    private FrequencyCalc fc;
    private double[][] cmi;
    private int nbAttrib;
    private int nbClasses;


    /**
     *  Constructor for the CondMutualInfo object
     *
     * @param  fc  Description of Parameter
     */
    public CondMutualInfo(FrequencyCalc fc) {
        this.names = fc.names;
        this.fc = fc;
        this.nbAttrib = names.length - 1;
        this.nbClasses = names[nbAttrib].getStates().length;
        this.cmi = new double[nbAttrib][nbAttrib];
    }


    /**
     *  Return conditional mutual information matrix I. The condition may be on
     *  the last variable, Z, in the sample data set. When both variables are
     *  conditioned on Z the following formula is used:
     * <pre>
     *              ___
     *             \                    p(x,y|z)
     *  I(X;Y|Z) = /__  p(x,y,z) log ---------------
     *             x,y                p(x|z) p(y|z)
     * </pre>
     * When only one variable, say X, is conditioned on Z:
     * <pre>
     *              ___
     *             \                   p(x,y|z)
     *  I(X;Y|Z) = /__  p(x,y,z) log -------------
     *             x,y                p(x|z) p(y)
     * </pre>
     * When neither variable is conditioned on Z: <pre>
     *            ___
     *           \                 p(x,y)
     *  I(X;Y) = /__  p(x,y) log -----------
     *           x,y              p(x) p(y)
     * </pre>
     *
     * @param  cond           weather a variable with corresponding index is
     *      conditioned on the class variable (Z).
     * @param  usePriors      Description of Parameter
     * @param  alphaK         Description of Parameter
     * @return                Description of the Returned Value
     * @exception  Exception  Description of Exception
     */
    public double[][] get(boolean[] cond,
                          boolean usePriors,
                          double alphaK) throws Exception {
        if (usePriors) {
            if (alphaK <= 0) {
                throw new Exception("alphaK has to be greater than zero.");
            }
        } else {
            alphaK = 0;
        }

        int[] n_z = fc.freqX[nbAttrib];

        for (int xi = 0; xi < nbAttrib; ++xi) {
            int size_x = names[xi].getStates().length;

            int[][][][] n_xYZ = fc.freqXYZ[xi];
            int[][][] n_xY = fc.freqXY[xi];
            int[][] n_xz = fc.freqXY[xi][nbAttrib];
            int[] n_x = fc.freqX[xi];

            double alpha_X = alphaK * size_x;

            cmi[xi][xi] = 0;

            for (int yi = xi + 1; yi < nbAttrib; ++yi) {
                int size_y = names[yi].getStates().length;

                int[][][] n_xyZ = n_xYZ[yi];
                int[][] n_xy = n_xY[yi];
                int[][] n_yz = fc.freqXY[yi][nbAttrib];
                int[] n_y = fc.freqX[yi];

                double alpha_Y = alphaK * size_y;
                double alpha_XY = alphaK * size_x * size_y;

                double v = 0;
                for (int vx = 0; vx < size_x; ++vx) {
                    int[][] n_xyZ_i = n_xyZ[vx];
                    int[] n_xy_i = n_xy[vx];
                    int[] n_xz_i = n_xz[vx];
                    int n_x_i = n_x[vx];
                    for (int vy = 0; vy < size_y; ++vy) {
                        int[] n_xyZ_ij = n_xyZ_i[vy];
                        int n_xy_ij = n_xy_i[vy];
                        int[] n_yz_j = n_yz[vy];
                        int n_y_j = n_y[vy];

                        if (!cond[xi] && !cond[yi]) {
                            // p(x), p(z)
                            if (!usePriors && (n_xy_ij == 0 || n_x_i == 0 || n_y_j == 0)) {
                                continue;
                            }

                            double a = (n_xy_ij + alphaK) * (fc.nbCases + alpha_X)
                                    * (fc.nbCases + alpha_Y);
                            a /= (fc.nbCases + alpha_XY) * (n_x_i + alphaK) * (n_y_j + alphaK);
                            double b = Math.log(a);
                            double m = (n_xy_ij + alphaK) / (fc.nbCases + alpha_XY) * b;
                            v += m;
                        } else {
                            double alpha_Z = alphaK * nbClasses;
                            double alpha_XZ = alphaK * size_x * nbClasses;
                            double alpha_YZ = alphaK * size_y * nbClasses;
                            double alpha_XYZ = alphaK * size_x * size_y * nbClasses;

                            for (int vz = 0; vz < nbClasses; ++vz) {
                                int n_xyZ_ijk = n_xyZ_ij[vz];
                                int n_xz_ik = n_xz_i[vz];
                                int n_yz_jk = n_yz_j[vz];
                                int n_z_k = n_z[vz];

                                if (n_xyZ_ijk == 0 || n_z_k == 0) {
                                    continue;
                                }

                                double a;

                                if (cond[xi] && cond[yi]) {
                                    // p(x|z), p(y|z)
                                    if (!usePriors && (n_xz_ik == 0 || n_yz_jk == 0)) {
                                        continue;
                                    }

                                    a = (n_xyZ_ijk + alphaK) * (n_z_k + alphaK)
                                            * (fc.nbCases + alpha_XZ) * (fc.nbCases + alpha_YZ);
                                    a /= (fc.nbCases + alpha_XYZ) * (fc.nbCases + alpha_Z)
                                            * (n_xz_ik + alphaK) * (n_yz_jk + alphaK);
                                } else if (cond[xi] && !cond[yi]) {
                                    // p(x|z), p(y)
                                    if (!usePriors && (n_xz_ik == 0 || n_y_j == 0)) {
                                        continue;
                                    }

//                  a = (double)(nbCases*n_xyZ_ijk)/(n_xz_ik*n_y_j);
                                    a = (n_xyZ_ijk + alphaK) * (fc.nbCases + alpha_XZ) * (fc.nbCases + alpha_Y);
                                    a /= (fc.nbCases + alpha_XYZ) * (n_xz_ik + alphaK) * (n_y_j + alphaK);
                                } else {
                                    // !cond[xi] && cond[yi]
                                    // p(x), p(y|z)
                                    if (!usePriors && (n_yz_jk == 0 || n_x_i == 0)) {
                                        continue;
                                    }

//                  a = (double)(nbCases*n_xyZ_ijk)/(n_x_i*n_yz_jk);
                                    a = (n_xyZ_ijk + alphaK) * (fc.nbCases + alpha_X) * (fc.nbCases + alpha_YZ);
                                    a /= (fc.nbCases + alpha_XYZ) * (n_x_i + alphaK) * (n_yz_jk + alphaK);
                                }

                                double b = Math.log(a);
//                double m = n_xyZ_ijk/(double)nbCases*b;
                                double m = (n_xyZ_ijk + alphaK) / (fc.nbCases + alpha_XYZ) * b;
                                v += m;
                            }
                            // for vz
                        }
                        // if
                    }
                    // for vy
                }
                // for vx

                cmi[xi][yi] = cmi[yi][xi] = v * inv_log2;
            }
            // for yi
        }
        // for xi

        // Dump
        /*
         *  java.io.FileOutputStream out  = new java.io.FileOutputStream("cmi_dump.txt");
         *  java.io.PrintStream      pOut = new java.io.PrintStream(out);
         *  for(int i=0; i<cmi.length; ++i)
         *  {
         *  for(int j=0; j<cmi[i].length; ++j)
         *  {
         *  pOut.println("["+i+","+j+"] = "+cmi[i][j]);
         *  }
         *  }
         */
        return cmi;
    }


    /**
     *  Return copy of mariginal probability of class variable, p(X).
     *
     * @return                The ProbX value
     * @exception  Exception  Not initialized.
     */
    public double[][] getProbX() throws Exception {
        double[][] probX = new double[fc.freqX.length][];
        for (int i = 0; i < probX.length; ++i) {
            int[] freqXi = fc.freqX[i];
            double[] probXi = new double[freqXi.length];
            for (int j = 0; j < probXi.length; ++j) {
                probXi[j] = (double) freqXi[j] / fc.nbCases;
            }

            probX[i] = probXi;
        }

        return probX;
    }
}