Changeset 8473 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess

Timestamp:

08/10/12 20:42:30 (12 years ago)

Author:

gkronber

Message:

#1902 worked on GPR: added line chart, made parameters of mean and covariance functions readable, removed target variable scaling, moved noise hyperparameter for likelihood function to the end of the parameter list, added methods to calculate the predicted variance, removed limits for scale of covariance functions and introduced exception handling to catch non-spd or singular cov matrixes, implemented rational quadratic covariance function, added unit test case from GBML book (however it does not work as the book seemingly uses a noise-less likelihood function)

Location:

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

Files:

• CovarianceConst.cs (modified) (2 diffs)
• CovarianceNoise.cs (modified) (2 diffs)
• CovariancePeriodic.cs (modified) (2 diffs)
• CovarianceRQiso.cs (added)
• CovarianceSEard.cs (modified) (2 diffs)
• CovarianceSEiso.cs (modified) (2 diffs)
• GaussianProcessModel.cs (modified) (10 diffs)
• GaussianProcessRegressionModelCreator.cs (modified) (2 diffs)
• GaussianProcessRegressionSolution.cs (modified) (2 diffs)
• MeanConst.cs (modified) (1 diff)
• MeanLinear.cs (modified) (1 diff)

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

@@ -32 +32 @@
     [Storable]
     private double sf2;
+    public double Scale { get { return sf2; } }
 
     [StorableConstructor]
@@ -56 +57 @@
 
     public void SetParameter(double[] hyp) {
-      this.sf2 = Math.Min(1E6, Math.Exp(2 * hyp[0])); // upper limit for scale
+      this.sf2 = Math.Exp(2 * hyp[0]);
     }
     public void SetData(double[,] x) {

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

@@ -32 +32 @@
     [Storable]
     private double sf2;
+    public double Scale { get { return sf2; } }
 
     [StorableConstructor]
@@ -56 +57 @@
 
     public void SetParameter(double[] hyp) {
-      this.sf2 = Math.Min(1E6, Math.Exp(2 * hyp[0])); // upper limit for scale
+      this.sf2 = Math.Exp(2 * hyp[0]);
     }
     public void SetData(double[,] x) {

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

@@ -35 +35 @@
     [Storable]
     private double sf2;
+    public double Scale { get { return sf2; } }
     [Storable]
     private double l;
+    public double Length { get { return l; } }
     [Storable]
     private double p;
+    public double Period { get { return p; } }
 
     private bool symmetric;
@@ -74 +77 @@
       this.p = Math.Exp(hyp[1]);
       this.sf2 = Math.Exp(2 * hyp[2]);
-
-      sf2 = Math.Min(10E6, sf2); // upper limit for the scale
+      // sf2 = Math.Min(10E6, sf2); // upper limit for the scale
 
       sd = null;

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

@@ -36 +36 @@
     [Storable]
     private double sf2;
+    public double Scale { get { return sf2; } }
+
     [Storable]
     private double[] l;
+    public double[] Length {
+      get {
+        if (l == null) return new double[0];
+        var copy = new double[l.Length];
+        Array.Copy(l, copy, copy.Length);
+        return copy;
+      }
+    }
 
     private double[,] sd;
@@ -76 +86 @@
       this.l = hyp.Take(hyp.Length - 1).Select(Math.Exp).ToArray();
       this.sf2 = Math.Exp(2 * hyp[hyp.Length - 1]);
-      sf2 = Math.Min(10E6, sf2); // upper limit for the scale
+      // sf2 = Math.Min(10E6, sf2); // upper limit for the scale
 
       sd = null;

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

@@ -37 +37 @@
     [Storable]
     private double sf2;
+    public double Scale { get { return sf2; } }
     [Storable]
     private double l;
+    public double Length { get { return l; } }
     [Storable]
     private bool symmetric;
@@ -80 +82 @@
     public void SetParameter(double[] hyp) {
       this.l = Math.Exp(hyp[0]);
-      this.sf2 = Math.Min(1E6, Math.Exp(2 * hyp[1])); // upper limit for scale
+      this.sf2 = Math.Exp(2 * hyp[1]);
       sd = null;
     }

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

@@ -74 +74 @@
     [Storable]
     private Scaling inputScaling;
-    [Storable]
-    private Scaling targetScaling;
 
 
@@ -85 +83 @@
       this.covarianceFunction = cloner.Clone(original.covarianceFunction);
       this.inputScaling = cloner.Clone(original.inputScaling);
-      this.targetScaling = cloner.Clone(original.targetScaling);
       this.negativeLogLikelihood = original.negativeLogLikelihood;
       this.targetVariable = original.targetVariable;
@@ -106 +103 @@
       this.allowedInputVariables = allowedInputVariables.ToArray();
 
-      sqrSigmaNoise = Math.Exp(2.0 * hyp.First());
-      sqrSigmaNoise = Math.Max(10E-6, sqrSigmaNoise); // lower limit for the noise level
 
       int nVariables = this.allowedInputVariables.Length;
-      this.meanFunction.SetParameter(hyp.Skip(1)
+      this.meanFunction.SetParameter(hyp
         .Take(this.meanFunction.GetNumberOfParameters(nVariables))
         .ToArray());
-      this.covarianceFunction.SetParameter(hyp.Skip(1 + this.meanFunction.GetNumberOfParameters(nVariables))
+      this.covarianceFunction.SetParameter(hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
         .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
         .ToArray());
+      sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
 
       CalculateModel(ds, rows);
@@ -123 +119 @@
       inputScaling = new Scaling(ds, allowedInputVariables, rows);
       x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
-
-
-      targetScaling = new Scaling(ds, new string[] { targetVariable }, rows);
-      var y = targetScaling.GetScaledValues(ds, targetVariable, rows);
+      var y = ds.GetDoubleValues(targetVariable, rows);
 
       int n = x.GetLength(0);
@@ -149 +142 @@
 
       var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
-      if (!res) throw new InvalidOperationException("Matrix is not positive semidefinite");
+      if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
 
       // calculate sum of diagonal elements for likelihood
@@ -198 +191 @@
       }
 
-      return new double[] { noiseGradient }
-        .Concat(meanGradients)
-        .Concat(covGradients).ToArray();
+      return
+        meanGradients
+        .Concat(covGradients)
+        .Concat(new double[] { noiseGradient }).ToArray();
     }
 
@@ -231 +225 @@
       // var kss = new double[newN];
       var Ks = new double[newN, n];
-      double[,] sWKs = new double[n, newN];
+      //double[,] sWKs = new double[n, newN];
       // double[,] v;
 
@@ -243 +237 @@
       var ms = meanFunction.GetMean(newX);
       for (int i = 0; i < newN; i++) {
-
         for (int j = 0; j < n; j++) {
           Ks[i, j] = covarianceFunction.GetCovariance(j, i);
-          sWKs[j, i] = Ks[i, j] / Math.Sqrt(sqrSigmaNoise);
+          //sWKs[j, i] = Ks[i, j] / Math.Sqrt(sqrSigmaNoise);
         }
       }
@@ -253 +246 @@
       // alglib.rmatrixsolvem(l, n, sWKs, newN, true, out info, out denseSolveRep, out v);
 
-      double targetScaleMin, targetScaleMax;
-      targetScaling.GetScalingParameters(targetVariable, out targetScaleMin, out targetScaleMax);
       return Enumerable.Range(0, newN)
-        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha))
-        .Select(m => m * (targetScaleMax - targetScaleMin) + targetScaleMin);
+        .Select(i => ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha));
       //for (int i = 0; i < newN; i++) {
       //  // predMean[i] = ms[i] + prod(GetRow(Ks, i), alpha);
@@ -265 +255 @@
 
     }
+
+    public IEnumerable<double> GetEstimatedVariance(Dataset dataset, IEnumerable<int> rows) {
+      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
+      int newN = newX.GetLength(0);
+      int n = x.GetLength(0);
+
+      var kss = new double[newN];
+      double[,] sWKs = new double[n, newN];
+
+
+      // for stddev 
+      covarianceFunction.SetData(newX);
+      for (int i = 0; i < newN; i++)
+        kss[i] = covarianceFunction.GetCovariance(i, i);
+
+      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 stddev 
+      int info;
+      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);
+
+      for (int i = 0; i < newN; i++) {
+        var sumV2 = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
+        yield return kss[i] - sumV2;
+      }
+    }
   }
 }

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

@@ -20 +20 @@
 #endregion
 
+using System;
 using System.Linq;
 using HeuristicLab.Common;
@@ -60 +61 @@
 
     public override IOperation Apply() {
-      var model = Create(ProblemData, Hyperparameter.ToArray(), MeanFunction, CovarianceFunction);
-      ModelParameter.ActualValue = model;
-      NegativeLogLikelihoodParameter.ActualValue = new DoubleValue(model.NegativeLogLikelihood);
-      HyperparameterGradientsParameter.ActualValue = new RealVector(model.GetHyperparameterGradients());
+      try {
+        var model = Create(ProblemData, Hyperparameter.ToArray(), MeanFunction, CovarianceFunction);
+        ModelParameter.ActualValue = model;
+        NegativeLogLikelihoodParameter.ActualValue = new DoubleValue(model.NegativeLogLikelihood);
+        HyperparameterGradientsParameter.ActualValue = new RealVector(model.GetHyperparameterGradients());
+        return base.Apply();
+      }
+      catch (ArgumentException) { }
+      catch (alglib.alglibexception) { }
+      NegativeLogLikelihoodParameter.ActualValue = new DoubleValue(1E300);
+      HyperparameterGradientsParameter.ActualValue = new RealVector(Hyperparameter.Count());
       return base.Apply();
     }

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

@@ -20 +20 @@
 #endregion
 
+using System.Collections.Generic;
+using System.Linq;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
@@ -51 +53 @@
       return new GaussianProcessRegressionSolution(this, cloner);
     }
+
+    public IEnumerable<double> EstimatedVariance {
+      get { return GetEstimatedVariance(Enumerable.Range(0, ProblemData.Dataset.Rows)); }
+    }
+    public IEnumerable<double> EstimatedTrainingVariance {
+      get { return GetEstimatedVariance(ProblemData.TrainingIndices); }
+    }
+    public IEnumerable<double> EstimatedTestVariance {
+      get { return GetEstimatedVariance(ProblemData.TestIndices); }
+    }
+
+    public IEnumerable<double> GetEstimatedVariance(IEnumerable<int> rows) {
+      return Model.GetEstimatedVariance(ProblemData.Dataset, rows);
+    }
   }
 }

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

@@ -32 +32 @@
     [Storable]
     private double c;
+    public double Value { get { return c; } }
+
     public int GetNumberOfParameters(int numberOfVariables) {
       return 1;

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

@@ -31 +31 @@
     [Storable]
     private double[] alpha;
+    public double[] Weights {
+      get {
+        if (alpha == null) return new double[0];
+        var copy = new double[alpha.Length];
+        Array.Copy(alpha, copy, copy.Length);
+        return copy;
+      }
+    }
     public int GetNumberOfParameters(int numberOfVariables) {
       return numberOfVariables;