Interface for multivariate implementations of the + * conditional mutual information.
+ * + *+ * Intended usage of the child classes: + *
Sets the single set of observations to compute the PDFs from. + * Cannot be called in conjunction with + * {@link #startAddObservations()}/{@link #addObservations(double[], double[], double[])} / + * {@link #finaliseAddObservations()}.
+ * + * @param var1 multivariate observations for variable 1 + * (first index is time, second is variable number) + * @param var2 multivariate observations for variable 2 + * (first index is time, second is variable number) + * @param cond multivariate observations for the conditional + * (first index is time, second is variable number) + * @throws Exception + */ + public void setObservations(double[][] var1, double[][] var2, + double[][] cond) throws Exception; + + /** + *Sets the single set of observations to compute the PDFs from. + * Cannot be called in conjunction with + * {@link #startAddObservations()}/{@link #addObservations(double[], double[])} / + * {@link #finaliseAddObservations()}.
+ * + * @param var1 multivariate observations for variable 1 + * (first index is time, second is variable number) + * @param var2 multivariate observations for variable 2 + * (first index is time, second is variable number) + * @param cond multivariate observations for the conditional + * (first index is time, second is variable number) + * @param var1Valid time series (with time indices the same as var1) + * indicating whether var1 at that point is valid. + * @param var2Valid time series (with time indices the same as var2) + * indicating whether var2 at that point is valid. + * @param condValid time series (with time indices the same as cond) + * indicating whether cond at that point is valid. + */ + public void setObservations(double[][] var1, double[][] var2, + double[][] cond, + boolean[] var1Valid, boolean[] var2Valid, + boolean[] condValid) throws Exception; + + /** + *Sets the single set of observations to compute the PDFs from. + * Cannot be called in conjunction with + * {@link #startAddObservations()}/{@link #addObservations(double[], double[])} / + * {@link #finaliseAddObservations()}.
+ * + * @param var1 multivariate observations for variable 1 + * (first index is time, second is variable number) + * @param var2 multivariate observations for variable 2 + * (first index is time, second is variable number) + * @param cond multivariate observations for the conditional + * (first index is time, second is variable number) + * @param var1Valid time series (with time indices the same as var1) + * indicating whether each variable of var1 at that point is valid. + * @param var2Valid time series (with time indices the same as var2) + * indicating whether each variable of var2 at that point is valid. + * @param condValid time series (with time indices the same as cond) + * indicating whether each variable of cond at that point is valid. + */ + public void setObservations(double[][] var1, double[][] var2, + double[][] cond, + boolean[][] var1Valid, boolean[][] var2Valid, + boolean[][] condValid) throws Exception; + + /** + * Elect to add in the observations from several disjoint time series. + * + */ + public void startAddObservations(); + + /** + *Adds a new set of observations to update the PDFs with - is + * intended to be called multiple times. + * Must be called after {@link #startAddObservations()}; call + * {@link #finaliseAddObservations()} once all observations have + * been supplied.
+ * + *Note that the arrays must not be over-written by the user + * until after finaliseAddObservations() has been called + * (they are not copied by this method necessarily, but the method + * may simply hold a pointer to them).
+ * + * @param var1 multivariate observations for variable 1 + * (first index is time, second is variable number) + * @param var2 multivariate observations for variable 2 + * (first index is time, second is variable number) + * @param cond multivariate observations for the conditional + * (first index is time, second is variable number) + * @throws Exception + */ + public void addObservations(double[][] var1, double[][] var2, + double[][] cond) throws Exception; + + /** + *Adds a new set of observations to update the PDFs with - is + * intended to be called multiple times. + * Must be called after {@link #startAddObservations()}; call + * {@link #finaliseAddObservations()} once all observations have + * been supplied.
+ * + *Note that the arrays must not be over-written by the user + * until after finaliseAddObservations() has been called + * (they are not copied by this method necessarily, but the method + * may simply hold a pointer to them).
+ * + * @param var1 multivariate observations for variable 1 + * (first index is time, second is variable number) + * @param var2 multivariate observations for variable 2 + * (first index is time, second is variable number) + * @param cond multivariate observations for the conditional + * (first index is time, second is variable number) + * @param startTime first time index to take observations on + * @param numTimeSteps number of time steps to use + * @throws Exception + */ + public void addObservations(double[][] var1, double[][] var2, + double[][] cond, + int startTime, int numTimeSteps) throws Exception; + + /** + * Flag that the observations are complete, probability distribution functions can now be built. + * @throws Exception + * + */ + public void finaliseAddObservations() throws Exception; + + /** + * + * @return the average value of the conditional mutual information measure, + * computed using all of the previously supplied observation sets. + * @throws Exception + */ + public double computeAverageLocalOfObservations() throws Exception; + + /** + *Computes the local values of the conditional mutual information, + * for each valid observation in the previously supplied observations + * (with PDFs computed using all of the previously supplied observation sets).
+ * + *If disjoint observations were supplied using several + * calls such as {@link ChannelCalculator#addObservations(double[], double[])} + * then the local values for each disjoint observation set will be appended here + * to create a single return array, + * though of course the time series for these disjoint observations were + * not appended in computing the required PDFs).
+ * + * @return array of local values. + * @throws Exception + */ + public double[] computeLocalOfPreviousObservations() throws Exception; + + /** + *Compute the significance of obtaining the given average + * measure from the given observations.
+ * + *This is in the spirit of Chavez et. al., "Statistical assessment of nonlinear causality: + * application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128. + *
+ * + *Basically, we shuffle the observations of the named variable + * against the other tuples. + * This keeps the marginal and joint PDFs of the unshuffled variables the same + * but destroys any correlation between the named variable and the others. + *
+ * + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param numPermutationsToCheck number of new orderings of the source values to compare against + * @see "Chavez et. al., 'Statistical assessment of nonlinear causality: + * application to epileptic EEG signals', Journal of Neuroscience Methods 124 (2003) 113-128" + */ + public EmpiricalMeasurementDistribution computeSignificance(int variableToReorder, + int numPermutationsToCheck) throws Exception; + + /** + *As per {@link #computeSignificance(int, int)} but supplies + * the re-orderings of the observations of the named variable.
+ * + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param newOrderings first index is permutation number, i.e. newOrderings[i] + * is an array of 1 permutation of 0..n-1, where there were n observations. + * If the length of each permutation in newOrderings + * is not equal to numObservations, an Exception is thrown. + * @return + * @throws Exception + */ + public EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int[][] newOrderings) throws Exception; + + /** + * Compute the conditional mutual information if the given variable were ordered as per the ordering + * specified in newOrdering + * + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param newOrdering permutation of the indices for the given variable + * @return + * @throws Exception + */ + public double computeAverageLocalOfObservations(int variableToReorder, int[] newOrdering) throws Exception; + + /** + * Compute the local mutual information for the given states, using the + * PDFs from the previously supplied observations. + * + * @param states1 + * @param states2 + * @param condStates + * @return + * @throws Exception + */ + public double[] computeLocalUsingPreviousObservations(double states1[][], double states2[][], double[][] condStates) + throws Exception; + + /** + * Set whether to print debug messages or not + * + * @param debug whether to print debug messages or not + */ + public void setDebug(boolean debug); + + /** + * Get the last computed average of the measure + * + * @return the last computed average + */ + public double getLastAverage(); + + /** + * Get the number of observations that have been supplied for + * computation of the PDFs + * + * @return number of observations + * @throws Exception + */ + public int getNumObservations() throws Exception; +} diff --git a/java/source/infodynamics/measures/continuous/ConditionalMutualInfoMultiVariateCommon.java b/java/source/infodynamics/measures/continuous/ConditionalMutualInfoMultiVariateCommon.java new file mode 100755 index 0000000..bfd60db --- /dev/null +++ b/java/source/infodynamics/measures/continuous/ConditionalMutualInfoMultiVariateCommon.java @@ -0,0 +1,484 @@ +package infodynamics.measures.continuous; + +import infodynamics.utils.EmpiricalMeasurementDistribution; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; + +import java.util.Iterator; +import java.util.Vector; + +/** + *Base class for implementations of {@link ConditionalMutualInfoCalculatorMultiVariate}, + * e.g. kernel estimation, Kraskov style extensions. + * It implements some common code to be used across conditional + * mutual information calculators + *
+ * + * + * @author Joseph Lizier, joseph.lizier at gmail.com + * + */ +public abstract class ConditionalMutualInfoMultiVariateCommon implements + ConditionalMutualInfoCalculatorMultiVariate { + + /** + * Number of dimenions for each of our multivariate data sets + */ + protected int dimensionsVar1 = 1; + protected int dimensionsVar2 = 1; + protected int dimensionsCond = 1; + + /** + * The set of observations for var1, retained in case the user wants to retrieve the local + * entropy values of these. + * They're held in the order in which they were supplied in the + * {@link #addObservations(double[][], double[][], double[][])} functions. + */ + protected double[][] var1Observations; + + /** + * The set of observations for var2, retained in case the user wants to retrieve the local + * entropy values of these. + * They're held in the order in which they were supplied in the + * {@link #addObservations(double[][], double[][], double[][])} functions. + */ + protected double[][] var2Observations; + + /** + * The set of observations for the conditional, retained in case the user wants to retrieve the local + * entropy values of these. + * They're held in the order in which they were supplied in the + * {@link #addObservations(double[][], double[][], double[][])} functions. + */ + protected double[][] condObservations; + + /** + * Total number of observations supplied. + * Only valid after {@link #finaliseAddObservations()} is called. + */ + protected int totalObservations = 0; + + /** + * Store the last computed average + */ + protected double lastAverage; + + /** + * Track whether we've computed the average for the supplied + * observations yet + */ + protected boolean condMiComputed; + + /** + * Whether to report debug messages or not + */ + protected boolean debug; + + /** + * Storage for var1 observations for addObservsations + */ + protected VectorAs per {@link #computeSignificance(int, int)} but supplies + * the re-orderings of the observations of the named variable.
+ * + *We provide a simple implementation which would be suitable for + * any child class, though the child class may prefer to make its + * own implementation to make class-specific optimisations. + * Child classes must implement {@link java.lang.Cloneable} + * for this method to be callable for them, and indeed implement + * the clone() method in a way that protects their structure + * from alteration by surrogate data being supplied to it.
+ * + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param newOrderings first index is permutation number, i.e. newOrderings[i] + * is an array of 1 permutation of 0..n-1, where there were n observations. + * If the length of each permutation in newOrderings + * is not equal to numObservations, an Exception is thrown. + * @param newOrderings the specific new orderings to use + * @return + * @throws Exception + */ + public EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int[][] newOrderings) throws Exception { + + int numPermutationsToCheck = newOrderings.length; + if (!condMiComputed) { + computeAverageLocalOfObservations(); + } + + // Take a clone of the object to compute the MI of the surrogates: + // (this is a shallow copy, it doesn't make new copies of all + // the arrays - child classes should override this) + ConditionalMutualInfoMultiVariateCommon miSurrogateCalculator = + (ConditionalMutualInfoMultiVariateCommon) this.clone(); + + double[] surrogateMeasurements = new double[numPermutationsToCheck]; + + // Now compute the MI for each set of shuffled data: + for (int i = 0; i < numPermutationsToCheck; i++) { + // Generate a new re-ordered source data + double[][] shuffledData = + MatrixUtils.extractSelectedTimePointsReusingArrays( + (variableToReorder == 1) ? var1Observations : var2Observations, + newOrderings[i]); + // Perform new initialisations + miSurrogateCalculator.initialise( + dimensionsVar1, dimensionsVar2, dimensionsCond); + // Set new observations + if (variableToReorder == 1) { + miSurrogateCalculator.setObservations(shuffledData, + var2Observations, condObservations); + } else { + miSurrogateCalculator.setObservations(var1Observations, + shuffledData, condObservations); + } + // Compute the MI + surrogateMeasurements[i] = miSurrogateCalculator.computeAverageLocalOfObservations(); + if (debug){ + System.out.println("New MI was " + surrogateMeasurements[i]); + } + } + + return new EmpiricalMeasurementDistribution(surrogateMeasurements, lastAverage); + } + + /** + *Compute the mutual information if the given variable were + * ordered as per the ordering specified in newOrdering.
+ * + *We provide a simple implementation which would be suitable for + * any child class, though the child class may prefer to make its + * own implementation to make class-specific optimisations. + * Child classes must implement {@link java.lang.Cloneable} + * for this method to be callable for them, and indeed implement + * the clone() method in a way that protects their structure + * from alteration by surrogate data being supplied to it.
+ * + * @param variableToReorder 1 for variable 1, 2 for variable 2 + * @param newOrdering array of time indices with which to reorder the data + * @return a surrogate MI evaluated for the given ordering of the source variable + * @throws Exception + */ + public double computeAverageLocalOfObservations(int variableToReorder, int[] newOrdering) + throws Exception { + // Take a clone of the object to compute the MI of the surrogates: + // (this is a shallow copy, it doesn't make new copies of all + // the arrays - child class should override this) + ConditionalMutualInfoMultiVariateCommon miSurrogateCalculator = + (ConditionalMutualInfoMultiVariateCommon) this.clone(); + + // Generate a new re-ordered source data + double[][] shuffledData = + MatrixUtils.extractSelectedTimePointsReusingArrays( + (variableToReorder == 1) ? var1Observations : var2Observations, + newOrdering); + // Perform new initialisations + miSurrogateCalculator.initialise( + dimensionsVar1, dimensionsVar2, dimensionsCond); + // Set new observations + if (variableToReorder == 1) { + miSurrogateCalculator.setObservations(shuffledData, + var2Observations, condObservations); + } else { + miSurrogateCalculator.setObservations(var1Observations, + shuffledData, condObservations); + } + // Compute the MI + return miSurrogateCalculator.computeAverageLocalOfObservations(); + } + + /** + * Set whether debug messages will be displayed + * + * @param debug debug setting + */ + public void setDebug(boolean debug) { + this.debug = debug; + } + + /** + * @return the previously computed average mutual information + */ + public double getLastAverage() { + return lastAverage; + } + + /** + * @return the number of supplied observations + * @throws Exception if child class computes MI without explicit observations + */ + public int getNumObservations() throws Exception { + return totalObservations; + } + + /** + * Compute a vector of start and end pairs of time points, between which we have + * valid series of all variables. + * + * Made public so it can be used if one wants to compute the number of + * observations prior to setting the observations. + * + * @param var1Valid + * @param var2Valid + * @return + */ + public VectorComputes the differential conditional mutual information of two given multivariate sets of + * observations, + * assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution.
+ * + *+ * Usage: + *
Set the covariance of the distribution for which we will compute the + * conditional mutual information.
+ * + *Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}, and without + * providing the means of the variables, you cannot later call + * {@link #computeLocalUsingPreviousObservations(double[][], double[][])}.
+ * + * @param covariance covariance matrix of var1, var2, conditional + * variables, considered together. + * @throws Exception for covariance matrix not matching the expected dimensions, + * being non-square, asymmetric or non-positive definite + */ + public void setCovariance(double[][] covariance) throws Exception { + setCovariance(covariance, false); + } + + /** + *Set the covariance of the distribution for which we will compute the + * conditional mutual information.
+ * + *Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}, and without + * providing the means of the variables, you cannot later call + * {@link #computeLocalUsingPreviousObservations(double[][], double[][])}.
+ * + * @param covariance covariance matrix of var1, var2 and the conditional + * variables, considered together. + * @param determinedFromObservations whether the covariance matrix + * was determined internally from observations or not + * @throws Exception for covariance matrix not matching the expected dimensions, + * being non-square, asymmetric or non-positive definite + */ + protected void setCovariance(double[][] covariance, boolean determinedFromObservations) + throws Exception { + if (!determinedFromObservations) { + // Make sure we're not keeping any observations + var1Observations = null; + var2Observations = null; + condObservations = null; + } + // Make sure the supplied covariance matrix matches the required dimenions: + int rows = covariance.length; + if (rows != dimensionsVar1 + dimensionsVar2 + dimensionsCond) { + throw new Exception("Supplied covariance matrix does not match initialised number of dimensions"); + } + + // Make sure the matrix is symmetric and positive definite, by taking the + // Cholesky decomposition (which we need for the determinant later anyway): + // (this will check and throw Exceptions for non-square, + // asymmetric, non-positive definite A) + L = MatrixUtils.CholeskyDecomposition(covariance); + + // Store the Cholesky decompositions for the conditional variable: + int[] condIndicesInCovariance = MatrixUtils.range(dimensionsVar1 + dimensionsVar2, + dimensionsVar1 + dimensionsVar2 + dimensionsCond - 1); + double[][] condCovariance = + MatrixUtils.selectRowsAndColumns(covariance, + condIndicesInCovariance, condIndicesInCovariance); + L_cc = MatrixUtils.CholeskyDecomposition(condCovariance); + // And store the Cholesky decompositions for var1 with + // the conditional variable: + int[] var1IndicesInCovariance = MatrixUtils.range(0, dimensionsVar1 - 1); + int[] var2IndicesInCovariance = MatrixUtils.range(dimensionsVar1, dimensionsVar1 + dimensionsVar2 - 1); + int[] var1AndCondIndicesInCovariance = MatrixUtils.append(var1IndicesInCovariance, condIndicesInCovariance); + int[] var2AndCondIndicesInCovariance = MatrixUtils.append(var2IndicesInCovariance, condIndicesInCovariance); + double[][] var1AndCondCovariance = + MatrixUtils.selectRowsAndColumns(covariance, + var1AndCondIndicesInCovariance, var1AndCondIndicesInCovariance); + L_1c = MatrixUtils.CholeskyDecomposition(var1AndCondCovariance); + double[][] var2AndCondCovariance = + MatrixUtils.selectRowsAndColumns(covariance, + var2AndCondIndicesInCovariance, var2AndCondIndicesInCovariance); + L_2c = MatrixUtils.CholeskyDecomposition(var2AndCondCovariance); + } + + /** + *Set the covariance of the distribution for which we will compute the + * mutual information.
+ * + *Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}.
+ * + * @param covariance covariance matrix of var1, var2 and conditional + * variables, considered together. + * @param means mean of var1, var2 and conditional variables (as per + * covariance) + */ + public void setCovarianceAndMeans(double[][] covariance, double[] means) throws Exception { + this.means = means; + setCovariance(covariance); + } + + /** + *The joint differential entropy for a multivariate Gaussian-distribution of dimension n + * with covariance matrix C is -0.5*\log_e{(2*pi*e)^n*|det(C)|}, + * where det() is the matrix determinant of C.
+ * + *Here we compute the conditional mutual information from the joint entropies + * of all variables (H_12c), variable 1 and conditional (H_1c), + * variable 2 and conditional (H_2c) and conditional (H_c), + * giving MI = H_1c + H_2c - H_c - H_12c. + * We assume that the recorded estimation of the + * covariance is correct (i.e. we will not make a bias correction for limited + * observations here).
+ * + * @return the mutual information of the previously provided observations or from the + * supplied covariance matrix, in nats (not bits!). + * Returns NaN if any of the determinants are zero + * (because this will make the denominator of the log zero) + */ + public double computeAverageLocalOfObservations() throws Exception { + // Simple way: + // detCovariance = MatrixUtils.determinantSymmPosDefMatrix(covariance); + // Using cached Cholesky decomposition: + detCovariance = MatrixUtils.determinantViaCholeskyResult(L); + det1cCovariance = MatrixUtils.determinantViaCholeskyResult(L_1c); + det2cCovariance = MatrixUtils.determinantViaCholeskyResult(L_2c); + detccCovariance = MatrixUtils.determinantViaCholeskyResult(L_cc); + + lastAverage = 0.5 * Math.log(Math.abs( + det1cCovariance * det2cCovariance / + (detCovariance * detccCovariance))); + condMiComputed = true; + return lastAverage; + } + + /** + *Compute the local or pointwise mutual information for each of the previously + * supplied observations
+ * + * @return array of the local values in nats (not bits!) + */ + public double[] computeLocalOfPreviousObservations() throws Exception { + // Cannot do if destObservations haven't been set + if (var2Observations == null) { + throw new Exception("Cannot compute local values of previous observations " + + "if they have not been set!"); + } + + return computeLocalUsingPreviousObservations(var1Observations, + var2Observations, condObservations, true); + } + + /** + *Compute the statistical significance of the conditional mutual information + * result analytically, without creating a distribution + * under the null hypothesis by bootstrapping.
+ * + *Brillinger (see reference below) shows that under the null hypothesis + * of no source-destination relationship, the MI for two + * Gaussian distributions follows a chi-square distribution with + * degrees of freedom equal to the product of the number of variables + * in each joint variable.
+ * + * @return ChiSquareMeasurementDistribution object + * This object contains the proportion of MI scores from the distribution + * which have higher or equal MIs to ours. + * + * @see Brillinger, "Some data analyses using mutual information", + * {@link http://www.stat.berkeley.edu/~brill/Papers/MIBJPS.pdf} + * @see Cheng et al., "Data Information in Contingency Tables: A + * Fallacy of Hierarchical Loglinear Models", + * {@link http://www.jds-online.com/file_download/112/JDS-369.pdf} + * @see Barnett and Bossomaier, "Transfer Entropy as a Log-likelihood Ratio" + * {@link http://arxiv.org/abs/1205.6339} + */ + public ChiSquareMeasurementDistribution computeSignificance() throws Exception { + if (!condMiComputed) { + computeAverageLocalOfObservations(); + } + // Number of extra parameters in the model incorporating the + // extra variable is independent of the number of variables + // in the conditional: + return new ChiSquareMeasurementDistribution(2*totalObservations*lastAverage, + dimensionsVar1 * dimensionsVar2); + } + + /** + * @return the number of previously supplied observations for which + * the conditional mutual information will be / was computed. + */ + public int getNumObservations() throws Exception { + if (var2Observations == null) { + throw new Exception("Cannot return number of observations because either " + + "this calculator has not had observations supplied or " + + "the user supplied the covariance matrix instead of observations"); + } + return super.getNumObservations(); + } + + /** + * Compute the conditional mutual information if the given variable was + * ordered as per the ordering specified in newOrdering + * + * @param newOrdering array of time indices with which to reorder the data + * @return a surrogate conditional MI evaluated for the given ordering of the source variable + * @throws Exception if the user previously supplied covariance directly rather + * than by setting observations (this means we have no observations + * to reorder). + */ + public double computeAverageLocalOfObservations(int variableToReorder, + int[] newOrdering) throws Exception { + // Cannot do if observations haven't been set (i.e. the variances + // were directly supplied) + if (var1Observations == null) { + throw new Exception("Cannot compute local values of previous observations " + + "without supplying observations"); + } + return super.computeAverageLocalOfObservations(variableToReorder, newOrdering); + } + + /** + * Compute the local conditional mutual information for a new series of + * observations, based on variances computed with the previously + * supplied observations. + * + * @param newVar1Obs provided variable 1 observations + * @param newVar2Obs provided variable 2 observations + * @param newCondObs provided conditional observations + * @return the local values in nats (not bits). + * @throws Exception + */ + public double[] computeLocalUsingPreviousObservations(double[][] newVar1Obs, + double[][] newVar2Obs, double[][] newCondObs) throws Exception { + return computeLocalUsingPreviousObservations( + newVar1Obs, newVar2Obs, newCondObs, false); + } + + /** + * Compute the local conditional mutual information for a new series of + * observations, based on variances computed with the previously + * supplied observations. + * + * @param newVar1Obs provided variable 1 observations + * @param newVar2Obs provided variable 2 observations + * @param newCondObs provided conditional observations + * @param isPreviousObservations whether these are our previous + * observations - this determines whether to + * set the internal lastAverage field, + * which is returned by later calls to {@link #getLastAverage()} + * @return the local values in nats (not bits). + * @see Multivariate normal distribution on Wikipedia + * @see + * @throws Exception if means were not defined by {@link #setObservations(double[][], double[][])} etc + * or {@link #setCovarianceAndMeans(double[][], double[])} + */ + protected double[] computeLocalUsingPreviousObservations(double[][] newVar1Obs, + double[][] newVar2Obs, double[][] newCondObs, boolean isPreviousObservations) throws Exception { + + if (means == null) { + throw new Exception("Cannot compute local values without having means either supplied or computed via setObservations()"); + } + + // Check that the covariance matrix was positive definite: + // (this was done earlier in computing the Cholesky decomposition, + // we may still need to compute the determinant) + if (detCovariance == 0) { + // The determinant has not been computed yet + // Simple way: + // detCovariance = MatrixUtils.determinantSymmPosDefMatrix(covariance); + // Using cached Cholesky decomposition: + detCovariance = MatrixUtils.determinantViaCholeskyResult(L); + if (detCovariance == 0) { + throw new Exception("Covariance matrix is not positive definite"); + } + det1cCovariance = MatrixUtils.determinantViaCholeskyResult(L_1c); + det2cCovariance = MatrixUtils.determinantViaCholeskyResult(L_2c); + detccCovariance = MatrixUtils.determinantViaCholeskyResult(L_cc); + } + + // Now we are clear to take the matrix inverse (via Cholesky decomposition, + // since we have a symmetric positive definite matrix): + double[][] invCovariance = MatrixUtils.solveViaCholeskyResult(L, + MatrixUtils.identityMatrix(L.length)); + double[][] invVar1CondCovariance = MatrixUtils.solveViaCholeskyResult(L_1c, + MatrixUtils.identityMatrix(L_1c.length)); + double[][] invVar2CondCovariance = MatrixUtils.solveViaCholeskyResult(L_2c, + MatrixUtils.identityMatrix(L_2c.length)); + double[][] invCondCovariance = MatrixUtils.solveViaCholeskyResult(L_cc, + MatrixUtils.identityMatrix(L_cc.length)); + + double[] var1Means = MatrixUtils.select(means, 0, dimensionsVar1); + double[] var2Means = MatrixUtils.select(means, dimensionsVar1, dimensionsVar2); + double[] condMeans = MatrixUtils.select(means, dimensionsVar1 + dimensionsVar2, dimensionsCond); + + int lengthOfReturnArray; + lengthOfReturnArray = newVar2Obs.length; + + double[] localValues = new double[lengthOfReturnArray]; + for (int t = 0; t < newVar2Obs.length; t++) { + + double[] var1DeviationsFromMean = + MatrixUtils.subtract(newVar1Obs[t], + var1Means); + double[] var2DeviationsFromMean = + MatrixUtils.subtract(newVar2Obs[t], var2Means); + double[] condDeviationsFromMean = + MatrixUtils.subtract(newCondObs[t], condMeans); + double[] var1CondDeviationsFromMean = + MatrixUtils.append(var1DeviationsFromMean, + condDeviationsFromMean); + double[] var2CondDeviationsFromMean = + MatrixUtils.append(var2DeviationsFromMean, + condDeviationsFromMean); + double[] tempDeviationsFromMean = + MatrixUtils.append(var1DeviationsFromMean, + var2DeviationsFromMean); + double[] deviationsFromMean = + MatrixUtils.append(tempDeviationsFromMean, + condDeviationsFromMean); + + // Computing PDFs WITHOUT (2*pi)^dim factor, since these will cancel: + // (see the PDFs defined at the wikipedia page referenced in the method header) + double var1CondExpArg = MatrixUtils.dotProduct( + MatrixUtils.matrixProduct(var1CondDeviationsFromMean, + invVar1CondCovariance), + var1CondDeviationsFromMean); + double adjustedPVar1Cond = Math.exp(-0.5 * var1CondExpArg) / + Math.sqrt(det1cCovariance); + double var2CondExpArg = MatrixUtils.dotProduct( + MatrixUtils.matrixProduct(var2CondDeviationsFromMean, + invVar2CondCovariance), + var2CondDeviationsFromMean); + double adjustedPVar2Cond = Math.exp(-0.5 * var2CondExpArg) / + Math.sqrt(det2cCovariance); + double condExpArg = MatrixUtils.dotProduct( + MatrixUtils.matrixProduct(condDeviationsFromMean, + invCondCovariance), + condDeviationsFromMean); + double adjustedPCond = Math.exp(-0.5 * condExpArg) / + Math.sqrt(detccCovariance); + double jointExpArg = MatrixUtils.dotProduct( + MatrixUtils.matrixProduct(deviationsFromMean, + invCovariance), + deviationsFromMean); + double adjustedPJoint = Math.exp(-0.5 * jointExpArg) / + Math.sqrt(detCovariance); + + // Returning results in nats: + double localValue = Math.log(adjustedPJoint * adjustedPCond / + (adjustedPVar1Cond * adjustedPVar2Cond)); + localValues[t] = localValue; + } + + // if (isPreviousObservations) { + // Don't store the average value here, since it won't be exactly + // the same as what would have been computed under the analytic expression + // } + + return localValues; + } + + /** + * No properties to set for this calculator + */ + public void setProperty(String propertyName, String propertyValue) + throws Exception { + + } + +}