Free cookie consent management tool by TermsFeed Policy Generator

Changeset 8475 for trunk


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis.Views/3.4/GaussianProcessRegressionSolutionLineChartView.cs

    r8473 r8475  
    6363        this.chart.ChartAreas[0].AxisX.Maximum = Content.ProblemData.Dataset.Rows - 1;
    6464
    65         this.chart.Series.Add(TARGETVARIABLE_SERIES_NAME);
    66         this.chart.Series[TARGETVARIABLE_SERIES_NAME].LegendText = Content.ProblemData.TargetVariable;
    67         this.chart.Series[TARGETVARIABLE_SERIES_NAME].ChartType = SeriesChartType.FastLine;
    68         this.chart.Series[TARGETVARIABLE_SERIES_NAME].Points.DataBindXY(Enumerable.Range(0, Content.ProblemData.Dataset.Rows).ToArray(),
    69           Content.ProblemData.Dataset.GetDoubleValues(Content.ProblemData.TargetVariable).ToArray());
    7065        // training series
    7166        this.chart.Series.Add(ESTIMATEDVALUES_TRAINING_SERIES_NAME);
     
    7570        var mean = Content.EstimatedTrainingValues.ToArray();
    7671        var s2 = Content.EstimatedTrainingVariance.ToArray();
    77         var lower = mean.Zip(s2, (m, s) => m - s).ToArray();
    78         var upper = mean.Zip(s2, (m, s) => m + s).ToArray();
     72        var lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
     73        var upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
    7974        this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME].Points.DataBindXY(Content.ProblemData.TrainingIndices.ToArray(), lower, upper);
    8075        this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME]);
    8176        this.chart.Series[ESTIMATEDVALUES_TRAINING_SERIES_NAME].Tag = Content;
     77
    8278        // test series
    8379        this.chart.Series.Add(ESTIMATEDVALUES_TEST_SERIES_NAME);
     
    8783        mean = Content.EstimatedTestValues.ToArray();
    8884        s2 = Content.EstimatedTestVariance.ToArray();
    89         lower = mean.Zip(s2, (m, s) => m - s).ToArray();
    90         upper = mean.Zip(s2, (m, s) => m + s).ToArray();
     85        lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
     86        upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
    9187        this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME].Points.DataBindXY(Content.ProblemData.TestIndices.ToArray(), lower, upper);
    9288        this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME]);
    9389        this.chart.Series[ESTIMATEDVALUES_TEST_SERIES_NAME].Tag = Content;
     90
    9491        // series of remaining points
    9592        int[] allIndices = Enumerable.Range(0, Content.ProblemData.Dataset.Rows).Except(Content.ProblemData.TrainingIndices).Except(Content.ProblemData.TestIndices).ToArray();
     
    10299        this.InsertEmptyPoints(this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME]);
    103100        this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME].Tag = Content;
     101
     102        // target
     103        this.chart.Series.Add(TARGETVARIABLE_SERIES_NAME);
     104        this.chart.Series[TARGETVARIABLE_SERIES_NAME].LegendText = Content.ProblemData.TargetVariable;
     105        this.chart.Series[TARGETVARIABLE_SERIES_NAME].ChartType = SeriesChartType.FastLine;
     106        this.chart.Series[TARGETVARIABLE_SERIES_NAME].Points.DataBindXY(Enumerable.Range(0, Content.ProblemData.Dataset.Rows).ToArray(),
     107          Content.ProblemData.Dataset.GetDoubleValues(Content.ProblemData.TargetVariable).ToArray());
     108
    104109        this.ToggleSeriesData(this.chart.Series[ESTIMATEDVALUES_ALL_SERIES_NAME]);
    105110
     
    243248            mean = Content.EstimatedValues.ToArray();
    244249            s2 = Content.EstimatedVariance.ToArray();
    245             lower = mean.Zip(s2, (m, s) => m - s).ToArray();
    246             upper = mean.Zip(s2, (m, s) => m + s).ToArray();
     250            lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
     251            upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
    247252            lower = indices.Select(index => lower[index]).ToArray();
    248253            upper = indices.Select(index => upper[index]).ToArray();
     
    252257            mean = Content.EstimatedTrainingValues.ToArray();
    253258            s2 = Content.EstimatedTrainingVariance.ToArray();
    254             lower = mean.Zip(s2, (m, s) => m - s).ToArray();
    255             upper = mean.Zip(s2, (m, s) => m + s).ToArray();
     259            lower = mean.Zip(s2, (m, s) => m - 1.96 * Math.Sqrt(s)).ToArray();
     260            upper = mean.Zip(s2, (m, s) => m + 1.96 * Math.Sqrt(s)).ToArray();
    256261            break;
    257262          case ESTIMATEDVALUES_TEST_SERIES_NAME:
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

    r8473 r8475  
    130130      double[] m = meanFunction.GetMean(x);
    131131      for (int i = 0; i < n; i++) {
    132 
    133132        for (int j = i; j < n; j++) {
    134133          l[j, i] = covarianceFunction.GetCovariance(i, j) / sqrSigmaNoise;
     
    163162      int info;
    164163      alglib.matinvreport matInvRep;
    165 
    166       alglib.spdmatrixcholeskyinverse(ref l, n, false, out info, out matInvRep);
     164      double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
     165      Array.Copy(l, lCopy, lCopy.Length);
     166
     167      alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
    167168      if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
    168169      for (int i = 0; i < n; i++) {
    169170        for (int j = 0; j <= i; j++)
    170           l[i, j] = l[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
    171       }
    172 
    173       double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => l[i, i]).Sum();
     171          lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
     172      }
     173
     174      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
    174175
    175176      double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
     
    184185          for (int k = 0; k < covGradients.Length; k++) {
    185186            for (int j = 0; j < i; j++) {
    186               covGradients[k] += l[i, j] * covarianceFunction.GetGradient(i, j, k);
     187              covGradients[k] += lCopy[i, j] * covarianceFunction.GetGradient(i, j, k);
    187188            }
    188             covGradients[k] += 0.5 * l[i, i] * covarianceFunction.GetGradient(i, i, k);
     189            covGradients[k] += 0.5 * lCopy[i, i] * covarianceFunction.GetGradient(i, i, k);
    189190          }
    190191        }
     
    264265      double[,] sWKs = new double[n, newN];
    265266
    266 
    267267      // for stddev
    268268      covarianceFunction.SetData(newX);
     
    271271
    272272      covarianceFunction.SetData(x, newX);
    273       for (int i = 0; i < n; i++) {
    274         for (int j = 0; j < newN; j++) {
    275           sWKs[i, j] = covarianceFunction.GetCovariance(i, j) / Math.Sqrt(sqrSigmaNoise);
     273      for (int i = 0; i < newN; i++) {
     274        for (int j = 0; j < n; j++) {
     275          sWKs[j, i] = covarianceFunction.GetCovariance(j, i) / Math.Sqrt(sqrSigmaNoise);
    276276        }
    277277      }
     
    281281      alglib.densesolverreport denseSolveRep;
    282282      double[,] v;
    283       double[,] lTrans = new double[l.GetLength(1), l.GetLength(0)];
    284       for (int i = 0; i < lTrans.GetLength(0); i++)
    285         for (int j = 0; j < lTrans.GetLength(1); j++)
    286           lTrans[i, j] = l[j, i];
    287       alglib.rmatrixsolvem(lTrans, n, sWKs, newN, true, out info, out denseSolveRep, out v); // not working!
    288       // alglib.spdmatrixcholeskysolvem(lTrans, n, true, sWKs, newN, out info, out denseSolveRep, out v);
     283
     284      alglib.rmatrixsolvem(l, n, sWKs, newN, false, out info, out denseSolveRep, out v);
    289285
    290286      for (int i = 0; i < newN; i++) {
    291         var sumV2 = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
    292         yield return kss[i] - sumV2;
    293       }
     287        var sumV = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
     288        kss[i] -= sumV;
     289        if (kss[i] < 0) kss[i] = 0;
     290      }
     291      return kss;
    294292    }
    295293  }
Note: See TracChangeset for help on using the changeset viewer.