Java tutorial
package compilateur.test_independance; import java.util.ArrayList; import org.apache.commons.math.stat.inference.ChiSquareTestImpl; import org.apache.commons.math.stat.inference.UnknownDistributionChiSquareTest; import compilateur.VDD; import compilateur.Var; /* (C) Copyright 2013, Schmidt Nicolas * * 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/>. */ public class TestG2 implements TestIndependance { protected UnknownDistributionChiSquareTest testStatistic = new ChiSquareTestImpl(); public double[][] getIndependancy(ArrayList<Var> v, VDD graph) { Var var1, var2; double[][] variance = new double[v.size()][v.size()]; for (int i = 0; i < v.size(); i++) { var1 = v.get(i); System.out.println(var1.name); for (int j = i; j < v.size(); j++) { //---debut du calcul----- var2 = v.get(j); variance[i][j] = computeInd(var1, var2, graph); System.out.print(var2.name + "=" + variance[i][j] + " "); } System.out.println(); } return variance; } private double computeInd(Var var1, Var var2, VDD graph) { int dom1, dom2; long[][] table; long[] sum1; long[] sum2; dom1 = var1.domain; dom2 = var2.domain; table = new long[dom1][dom2]; sum1 = new long[dom1]; sum2 = new long[dom2]; for (int l = 0; l < dom1; l++) sum1[l] = 0; for (int l = 0; l < dom2; l++) sum2[l] = 0; int N = 0; //calcul des proba au cas par cas for (int l = 0; l < dom1; l++) { graph.conditioner(var1, l); for (int k = 0; k < dom2; k++) { graph.conditioner(var2, k); table[l][k] = graph.countingpondere(); graph.deconditioner(var2); sum1[l] += table[l][k]; sum2[k] += table[l][k]; N += table[l][k]; // System.out.println("table["+l+"]["+k+"] : "+table[l][k]); } graph.deconditioner(var1); } int tailleV1 = dom1; int tailleV2 = dom2; double[] sommeV1 = new double[tailleV1]; double[] sommeV2 = new double[tailleV2]; boolean[] suppr1 = new boolean[tailleV1]; boolean[] suppr2 = new boolean[tailleV2]; int newTailleV1 = tailleV1; int newTailleV2 = tailleV2; for (int i = 0; i < tailleV1; i++) { int somme = 0; for (int j = 0; j < tailleV2; j++) somme += table[i][j]; sommeV1[i] = somme; if (somme == 0) { suppr1[i] = true; newTailleV1--; } else suppr1[i] = false; } for (int i = 0; i < tailleV2; i++) { int somme = 0; for (int j = 0; j < tailleV1; j++) somme += table[j][i]; sommeV2[i] = somme; if (somme == 0) { suppr2[i] = true; newTailleV2--; } else suppr2[i] = false; } long[][] newTable = new long[newTailleV1][newTailleV2]; long[] newSommeV1 = new long[newTailleV1]; long[] newSommeV2 = new long[newTailleV2]; int a = 0, b = 0; for (int i = 0; i < tailleV1; i++) { if (sommeV1[i] == 0) continue; int somme = 0; for (int j = 0; j < tailleV2; j++) somme += table[i][j]; newSommeV1[a] = somme; a++; } for (int i = 0; i < tailleV2; i++) { if (sommeV2[i] == 0) continue; int somme = 0; for (int j = 0; j < tailleV1; j++) somme += table[j][i]; newSommeV2[b] = somme; b++; } a = 0; b = 0; for (int i = 0; i < tailleV1; i++) if (suppr1[i]) continue; else { b = 0; for (int j = 0; j < tailleV2; j++) { if (suppr2[j]) continue; else newTable[a][b++] = table[i][j]; } a++; } // for(int l = 0; l < newTailleV1; l++) // for(int k = 0; k < newTailleV2; k++) // System.out.println("newTable["+l+"]["+k+"] : "+newTable[l][k]); /* try { return testStatistic.chiSquareTest(newTable); } catch (IllegalArgumentException e) { e.printStackTrace(); } catch (MathException e) { e.printStackTrace(); } return -1;*/ double statistique = 0; for (int l = 0; l < newTailleV1; l++) for (int k = 0; k < newTailleV2; k++) if (newTable[l][k] != 0) { statistique += 2 * newTable[l][k] * Math.log(((double) newTable[l][k] * N) / (newSommeV1[l] * newSommeV2[k])); /* System.out.println("N:"+N); System.out.println("newSommeV1:"+newSommeV1[l]); System.out.println("newSommeV2:"+newSommeV2[k]); System.out.println("newTable[l][k]:"+newTable[l][k]); System.out.println("statistique:"+statistique);*/ } // return statistique; return pochisq(statistique, (newTailleV1 - 1) * (newTailleV2 - 1)); // return g2(dom1, dom2, table, (dom1-1)*(dom2-1)); } private static final double LOG_SQRT_PI = Math.log(Math.sqrt(Math.PI)); private static final double I_SQRT_PI = 1 / Math.sqrt(Math.PI); private static final int MAX_X = 20; // max value to represent exp(x) /* POCHISQ -- probability of chi-square value Adapted from: Hill, I. D. and Pike, M. C. Algorithm 299 Collected Algorithms for the CACM 1967 p. 243 Updated for rounding errors based on remark in ACM TOMS June 1985, page 185 */ private double pochisq(double x, int df) { double a, s; double e, c, z; if (x <= 0.0 || df < 1) { return 1.0; } a = 0.5 * x; boolean even = (df & 1) == 0; double y = 0; if (df > 1) { y = Math.exp(-a); } s = (even ? y : (2.0 * poz(-Math.sqrt(x)))); if (df > 2) { x = 0.5 * (df - 1.0); z = (even ? 1.0 : 0.5); if (a > MAX_X) { e = (even ? 0.0 : LOG_SQRT_PI); c = Math.log(a); while (z <= x) { e = Math.log(z) + e; s += Math.exp(c * z - a - e); z += 1.0; } return s; } else { e = (even ? 1.0 : (I_SQRT_PI / Math.sqrt(a))); c = 0.0; while (z <= x) { e = e * (a / z); c = c + e; z += 1.0; } return c * y + s; } } else { return s; } } /* private double g2(int tailleV1, int tailleV2, double[][] table, int dfcorr) { double statistique = 0; double[] sommeV1 = new double[tailleV1]; double[] sommeV2 = new double[tailleV2]; boolean[] suppr1 = new boolean[tailleV1]; boolean[] suppr2 = new boolean[tailleV2]; int newTailleV1 = tailleV1; int newTailleV2 = tailleV2; for(int i = 0; i < tailleV1; i++) { int somme = 0; for(int j = 0; j < tailleV2; j++) somme += table[i][j]; sommeV1[i] = somme; if(somme == 0) { suppr1[i] = true; newTailleV1--; } else suppr1[i] = false; } for(int i = 0; i < tailleV2; i++) { int somme = 0; for(int j = 0; j < tailleV1; j++) somme += table[j][i]; sommeV2[i] = somme; if(somme == 0) { suppr2[i] = true; newTailleV2--; } else suppr2[i] = false; } double[][] newTable = new double[newTailleV1][newTailleV2]; double[] newSommeV1 = new double[newTailleV1]; double[] newSommeV2 = new double[newTailleV2]; int a = 0, b = 0; for(int i = 0; i < tailleV1; i++) { if(sommeV1[i] == 0) continue; int somme = 0; for(int j = 0; j < tailleV2; j++) somme += table[i][j]; newSommeV1[a] = somme; a++; } for(int i = 0; i < tailleV2; i++) { if(sommeV2[i] == 0) continue; int somme = 0; for(int j = 0; j < tailleV1; j++) somme += table[j][i]; newSommeV2[b] = somme; b++; } a = 0; b = 0; for(int i = 0; i < tailleV1; i++) if(suppr1[i]) continue; else { b = 0; for(int j = 0; j < tailleV2; j++) { if(suppr2[j]) continue; else newTable[a][b++] = table[i][j]; } a++; } int N = 0; for(int i = 0; i < tailleV1; i++) for(int j = 0; j < tailleV2; j++) N += table[i][j]; // pour la correction du degr de libert, voir Steck & Jaakkola (2002) for(int l=0; l<newTailleV1; l++) for(int k=0; k<newTailleV2; k++) if(table[l][k] != 0) statistique += 2*table[l][k]*Math.log(table[l][k]*N/(sommeV1[l]*sommeV2[k])); // return statistique; return pochisq(statistique, (newTailleV1-1)*(newTailleV2-1)*dfcorr); }*/ private double poz(double z) { double y, x, w; double Z_MAX = 6.0; // Maximum meaningful z value if (z == 0.0) { x = 0.0; } else { y = 0.5 * Math.abs(z); if (y >= (Z_MAX * 0.5)) { x = 1.0; } else if (y < 1.0) { w = y * y; x = ((((((((0.000124818987 * w - 0.001075204047) * w + 0.005198775019) * w - 0.019198292004) * w + 0.059054035642) * w - 0.151968751364) * w + 0.319152932694) * w - 0.531923007300) * w + 0.797884560593) * y * 2.0; } else { y -= 2.0; x = (((((((((((((-0.000045255659 * y + 0.000152529290) * y - 0.000019538132) * y - 0.000676904986) * y + 0.001390604284) * y - 0.000794620820) * y - 0.002034254874) * y + 0.006549791214) * y - 0.010557625006) * y + 0.011630447319) * y - 0.009279453341) * y + 0.005353579108) * y - 0.002141268741) * y + 0.000535310849) * y + 0.999936657524; } } return z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5); } @Override public boolean estPlusIndependantQue(double valeur1, double valeur2) { return valeur1 > valeur2; } public boolean isTestSymetrique() { return true; } }