Changeset 8475

Timestamp:

08/11/12 14:45:15 (12 years ago)

Author:

gkronber

Message:

#1902 fixed bug in calculation of variance in GPR model

Location:

trunk/sources

Files:

• HeuristicLab.Algorithms.DataAnalysis.Views/3.4/GaussianProcessRegressionSolutionLineChartView.cs (modified) (6 diffs)
• HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs (modified) (6 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis.Views/3.4/GaussianProcessRegressionSolutionLineChartView.cs

@@ -63 +63 @@
         this.chart.ChartAreas[0].AxisX.Maximum = Content.ProblemData.Dataset.Rows - 1;
 
-        this.chart.Series.Add(TARGETVARIABLE_SERIES_NAME);
-        this.chart.Series[TARGETVARIABLE_SERIES_NAME].LegendText = Content.ProblemData.TargetVariable;
-        this.chart.Series[TARGETVARIABLE_SERIES_NAME].ChartType = SeriesChartType.FastLine;
-        this.chart.Series[TARGETVARIABLE_SERIES_NAME].Points.DataBindXY(Enumerable.Range(0, Content.ProblemData.Dataset.Rows).ToArray(),
-          Content.ProblemData.Dataset.GetDoubleValues(Content.ProblemData.TargetVariable).ToArray());
         // training series
         this.chart.Series.Add(ESTIMATEDVALUES_TRAINING_SERIES_NAME);
@@ -75 +70 @@
         var mean = Content.EstimatedTrainingValues.ToArray();
         var s2 = Content.EstimatedTrainingVariance.ToArray();
-        var lower = mean.Zip(s2, (m, s) => m - s).ToArray();
-        var upper = mean.Zip(s2, (m, s) => m + s).ToArray();
+        var lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
+        var upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
         this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME].Points.DataBindXY(Content.ProblemData.TrainingIndices.ToArray(), lower, upper);
         this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME]);
         this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME].Tag = Content;
+
         // test series
         this.chart.Series.Add(ESTIMATEDVALUES_TEST_SERIES_NAME);
@@ -87 +83 @@
         mean = Content.EstimatedTestValues.ToArray();
         s2 = Content.EstimatedTestVariance.ToArray();
-        lower = mean.Zip(s2, (m, s) => m - s).ToArray();
-        upper = mean.Zip(s2, (m, s) => m + s).ToArray();
+        lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
+        upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
         this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME].Points.DataBindXY(Content.ProblemData.TestIndices.ToArray(), lower, upper);
         this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME]);
         this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME].Tag = Content;
+
         // series of remaining points
         int[] allIndices = Enumerable.Range(0, Content.ProblemData.Dataset.Rows).Except(Content.ProblemData.TrainingIndices).Except(Content.ProblemData.TestIndices).ToArray();
@@ -102 +99 @@
         this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME]);
         this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME].Tag = Content;
+
+        // target
+        this.chart.Series.Add(TARGETVARIABLE_SERIES_NAME);
+        this.chart.Series[TARGETVARIABLE_SERIES_NAME].LegendText = Content.ProblemData.TargetVariable;
+        this.chart.Series[TARGETVARIABLE_SERIES_NAME].ChartType = SeriesChartType.FastLine;
+        this.chart.Series[TARGETVARIABLE_SERIES_NAME].Points.DataBindXY(Enumerable.Range(0, Content.ProblemData.Dataset.Rows).ToArray(),
+          Content.ProblemData.Dataset.GetDoubleValues(Content.ProblemData.TargetVariable).ToArray());
+
         this.ToggleSeriesData(this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME]);
 
@@ -243 +248 @@
             mean = Content.EstimatedValues.ToArray();
             s2 = Content.EstimatedVariance.ToArray();
-            lower = mean.Zip(s2, (m, s) => m - s).ToArray();
-            upper = mean.Zip(s2, (m, s) => m + s).ToArray();
+            lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
+            upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
             lower = indices.Select(index => lower[index]).ToArray();
             upper = indices.Select(index => upper[index]).ToArray();
@@ -252 +257 @@
             mean = Content.EstimatedTrainingValues.ToArray();
             s2 = Content.EstimatedTrainingVariance.ToArray();
-            lower = mean.Zip(s2, (m, s) => m - s).ToArray();
-            upper = mean.Zip(s2, (m, s) => m + s).ToArray();
+            lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
+            upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
             break;
           case ESTIMATEDVALUES_TEST_SERIES_NAME:

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

@@ -130 +130 @@
       double[] m = meanFunction.GetMean(x);
       for (int i = 0; i < n; i++) {
-
         for (int j = i; j < n; j++) {
           l[j, i] = covarianceFunction.GetCovariance(i, j) / sqrSigmaNoise;
@@ -163 +162 @@
       int info;
       alglib.matinvreport matInvRep;
-
-      alglib.spdmatrixcholeskyinverse(ref l, n, false, out info, out matInvRep);
+      double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
+      Array.Copy(l, lCopy, lCopy.Length);
+
+      alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
       if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
       for (int i = 0; i < n; i++) {
         for (int j = 0; j <= i; j++)
-          l[i, j] = l[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
-      }
-
-      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => l[i, i]).Sum();
+          lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
+      }
+
+      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
 
       double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
@@ -184 +185 @@
           for (int k = 0; k < covGradients.Length; k++) {
             for (int j = 0; j < i; j++) {
-              covGradients[k] += l[i, j] * covarianceFunction.GetGradient(i, j, k);
+              covGradients[k] += lCopy[i, j] * covarianceFunction.GetGradient(i, j, k);
             }
-            covGradients[k] += 0.5 * l[i, i] * covarianceFunction.GetGradient(i, i, k);
+            covGradients[k] += 0.5 * lCopy[i, i] * covarianceFunction.GetGradient(i, i, k);
           }
         }
@@ -264 +265 @@
       double[,] sWKs = new double[n, newN];
 
-
       // for stddev 
       covarianceFunction.SetData(newX);
@@ -271 +271 @@
 
       covarianceFunction.SetData(x, newX);
-      for (int i = 0; i < n; i++) {
-        for (int j = 0; j < newN; j++) {
-          sWKs[i, j] = covarianceFunction.GetCovariance(i, j) / Math.Sqrt(sqrSigmaNoise);
+      for (int i = 0; i < newN; i++) {
+        for (int j = 0; j < n; j++) {
+          sWKs[j, i] = covarianceFunction.GetCovariance(j, i) / Math.Sqrt(sqrSigmaNoise);
         }
       }
@@ -281 +281 @@
       alglib.densesolverreport denseSolveRep;
       double[,] v;
-      double[,] lTrans = new double[l.GetLength(1), l.GetLength(0)];
-      for (int i = 0; i < lTrans.GetLength(0); i++)
-        for (int j = 0; j < lTrans.GetLength(1); j++)
-          lTrans[i, j] = l[j, i];
-      alglib.rmatrixsolvem(lTrans, n, sWKs, newN, true, out info, out denseSolveRep, out v); // not working!
-      // alglib.spdmatrixcholeskysolvem(lTrans, n, true, sWKs, newN, out info, out denseSolveRep, out v);
+
+      alglib.rmatrixsolvem(l, n, sWKs, newN, false, out info, out denseSolveRep, out v);
 
       for (int i = 0; i < newN; i++) {
-        var sumV2 = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
-        yield return kss[i] - sumV2;
-      }
+        var sumV = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
+        kss[i] -= sumV;
+        if (kss[i] < 0) kss[i] = 0;
+      }
+      return kss;
     }
   }