Changeset 8477 for branches/HeuristicLab.TimeSeries/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

Timestamp:

08/13/12 16:18:37 (12 years ago)

Author:

mkommend

Message:

#1081:

• Added autoregressive target variable Symbol
• Merged trunk changes into the branch.

Location:

branches/HeuristicLab.TimeSeries/HeuristicLab.Algorithms.DataAnalysis

Files:

• . (modified) (1 prop)
• 3.4/GaussianProcess/GaussianProcessModel.cs (modified) (12 diffs)

branches/HeuristicLab.TimeSeries/HeuristicLab.Algorithms.DataAnalysis/3.4/GaussianProcess/GaussianProcessModel.cs

@@ -73 +73 @@
     private double[,] x;
     [Storable]
-    private Scaling scaling;
+    private Scaling inputScaling;
 
 
@@ -82 +82 @@
       this.meanFunction = cloner.Clone(original.meanFunction);
       this.covarianceFunction = cloner.Clone(original.covarianceFunction);
-      this.scaling = cloner.Clone(original.scaling);
+      this.inputScaling = cloner.Clone(original.inputScaling);
       this.negativeLogLikelihood = original.negativeLogLikelihood;
       this.targetVariable = original.targetVariable;
@@ -103 +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);
@@ -118 +117 @@
 
     private void CalculateModel(Dataset ds, IEnumerable<int> rows) {
-      scaling = new Scaling(ds, allowedInputVariables, rows);
-      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, scaling);
-
-      var y = ds.GetDoubleValues(targetVariable, rows).ToArray();
+      inputScaling = new Scaling(ds, allowedInputVariables, rows);
+      x = AlglibUtil.PrepareAndScaleInputMatrix(ds, allowedInputVariables, rows, inputScaling);
+      var y = ds.GetDoubleValues(targetVariable, rows);
 
       int n = x.GetLength(0);
@@ -132 +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;
@@ -144 +141 @@
 
       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
@@ -162 +159 @@
       int n = x.GetLength(0);
       int nAllowedVariables = x.GetLength(1);
-      double[,] q = new double[n, n];
-      double[,] eye = new double[n, n];
-      for (int i = 0; i < n; i++) eye[i, i] = 1.0;
 
       int info;
-      alglib.densesolverreport denseSolveRep;
-
-      alglib.spdmatrixcholeskysolvem(l, n, false, eye, n, out info, out denseSolveRep, out q);
-      // double[,] a2 = outerProd(alpha, alpha);
+      alglib.matinvreport 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 < n; j++)
-          q[i, j] = q[i, j] / sqrSigmaNoise - alpha[i] * alpha[j]; // a2[i, j];
-      }
-
-      double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => q[i, i]).Sum();
+        for (int j = 0; j <= i; j++)
+          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)];
@@ -187 +183 @@
       if (covGradients.Length > 0) {
         for (int i = 0; i < n; i++) {
-          for (int j = 0; j < n; j++) {
-            var covDeriv = covarianceFunction.GetGradient(i, j);
-            for (int k = 0; k < covGradients.Length; k++) {
-              covGradients[k] += q[i, j] * covDeriv[k];
+          for (int k = 0; k < covGradients.Length; k++) {
+            for (int j = 0; j < i; j++) {
+              covGradients[k] += lCopy[i, j] * covarianceFunction.GetGradient(i, j, k);
             }
+            covGradients[k] += 0.5 * lCopy[i, i] * covarianceFunction.GetGradient(i, i, k);
           }
         }
-        covGradients = covGradients.Select(g => g / 2.0).ToArray();
-      }
-
-      return new double[] { noiseGradient }
-        .Concat(meanGradients)
-        .Concat(covGradients).ToArray();
+      }
+
+      return
+        meanGradients
+        .Concat(covGradients)
+        .Concat(new double[] { noiseGradient }).ToArray();
     }
 
@@ -220 +216 @@
 
     private IEnumerable<double> GetEstimatedValuesHelper(Dataset dataset, IEnumerable<int> rows) {
-      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, scaling);
+      var newX = AlglibUtil.PrepareAndScaleInputMatrix(dataset, allowedInputVariables, rows, inputScaling);
       int newN = newX.GetLength(0);
       int n = x.GetLength(0);
@@ -230 +226 @@
       // 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;
 
@@ -242 +238 @@
       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);
         }
       }
@@ -252 +247 @@
       // alglib.rmatrixsolvem(l, n, sWKs, newN, true, out info, out denseSolveRep, out v);
 
-
+      return Enumerable.Range(0, newN)
+        .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);
+      //  // var sumV2 = prod(GetCol(v, i), GetCol(v, i));
+      //  // predVar[i] = kss[i] - sumV2;
+      //}
+
+    }
+
+    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 < newN; i++) {
-        // predMean[i] = ms[i] + prod(GetRow(Ks, i), alpha);
-        yield return ms[i] + Util.ScalarProd(Util.GetRow(Ks, i), alpha);
-        // var sumV2 = prod(GetCol(v, i), GetCol(v, i));
-        // predVar[i] = kss[i] - sumV2;
-      }
-
+        for (int j = 0; j < n; j++) {
+          sWKs[j, i] = covarianceFunction.GetCovariance(j, i) / Math.Sqrt(sqrSigmaNoise);
+        }
+      }
+
+      // for stddev 
+      int info;
+      alglib.densesolverreport denseSolveRep;
+      double[,] v;
+
+      alglib.rmatrixsolvem(l, n, sWKs, newN, false, out info, out denseSolveRep, out v);
+
+      for (int i = 0; i < newN; i++) {
+        var sumV = Util.ScalarProd(Util.GetCol(v, i), Util.GetCol(v, i));
+        kss[i] -= sumV;
+        if (kss[i] < 0) kss[i] = 0;
+      }
+      return kss;
     }
   }