[8323] | 1 | #region License Information
|
---|
| 2 | /* HeuristicLab
|
---|
[14185] | 3 | * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
|
---|
[8323] | 4 | *
|
---|
| 5 | * This file is part of HeuristicLab.
|
---|
| 6 | *
|
---|
| 7 | * HeuristicLab is free software: you can redistribute it and/or modify
|
---|
| 8 | * it under the terms of the GNU General Public License as published by
|
---|
| 9 | * the Free Software Foundation, either version 3 of the License, or
|
---|
| 10 | * (at your option) any later version.
|
---|
| 11 | *
|
---|
| 12 | * HeuristicLab is distributed in the hope that it will be useful,
|
---|
| 13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 15 | * GNU General Public License for more details.
|
---|
| 16 | *
|
---|
| 17 | * You should have received a copy of the GNU General Public License
|
---|
| 18 | * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 19 | */
|
---|
| 20 | #endregion
|
---|
| 21 |
|
---|
| 22 | using System;
|
---|
| 23 | using System.Collections.Generic;
|
---|
| 24 | using System.Linq;
|
---|
| 25 | using HeuristicLab.Common;
|
---|
| 26 | using HeuristicLab.Core;
|
---|
[14927] | 27 | using HeuristicLab.Persistence;
|
---|
[8323] | 28 | using HeuristicLab.Problems.DataAnalysis;
|
---|
| 29 |
|
---|
[8371] | 30 | namespace HeuristicLab.Algorithms.DataAnalysis {
|
---|
[8323] | 31 | /// <summary>
|
---|
| 32 | /// Represents a Gaussian process model.
|
---|
| 33 | /// </summary>
|
---|
[14927] | 34 | [StorableType("36ca62fa-4766-4269-b5f1-3bc24c0ed2e1")]
|
---|
[8323] | 35 | [Item("GaussianProcessModel", "Represents a Gaussian process posterior.")]
|
---|
[13941] | 36 | public sealed class GaussianProcessModel : RegressionModel, IGaussianProcessModel {
|
---|
| 37 | public override IEnumerable<string> VariablesUsedForPrediction {
|
---|
[13922] | 38 | get { return allowedInputVariables; }
|
---|
| 39 | }
|
---|
[13921] | 40 |
|
---|
[8323] | 41 | [Storable]
|
---|
| 42 | private double negativeLogLikelihood;
|
---|
| 43 | public double NegativeLogLikelihood {
|
---|
| 44 | get { return negativeLogLikelihood; }
|
---|
| 45 | }
|
---|
| 46 |
|
---|
| 47 | [Storable]
|
---|
[14899] | 48 | private double negativeLooPredictiveProbability;
|
---|
| 49 | public double NegativeLooPredictiveProbability {
|
---|
| 50 | get { return negativeLooPredictiveProbability; }
|
---|
| 51 | }
|
---|
| 52 |
|
---|
| 53 | [Storable]
|
---|
[8484] | 54 | private double[] hyperparameterGradients;
|
---|
| 55 | public double[] HyperparameterGradients {
|
---|
| 56 | get {
|
---|
| 57 | var copy = new double[hyperparameterGradients.Length];
|
---|
| 58 | Array.Copy(hyperparameterGradients, copy, copy.Length);
|
---|
| 59 | return copy;
|
---|
| 60 | }
|
---|
| 61 | }
|
---|
| 62 |
|
---|
| 63 | [Storable]
|
---|
[8323] | 64 | private ICovarianceFunction covarianceFunction;
|
---|
| 65 | public ICovarianceFunction CovarianceFunction {
|
---|
| 66 | get { return covarianceFunction; }
|
---|
| 67 | }
|
---|
| 68 | [Storable]
|
---|
| 69 | private IMeanFunction meanFunction;
|
---|
| 70 | public IMeanFunction MeanFunction {
|
---|
| 71 | get { return meanFunction; }
|
---|
| 72 | }
|
---|
[13941] | 73 |
|
---|
[8323] | 74 | [Storable]
|
---|
| 75 | private string[] allowedInputVariables;
|
---|
| 76 | public string[] AllowedInputVariables {
|
---|
| 77 | get { return allowedInputVariables; }
|
---|
| 78 | }
|
---|
| 79 |
|
---|
| 80 | [Storable]
|
---|
| 81 | private double[] alpha;
|
---|
| 82 | [Storable]
|
---|
| 83 | private double sqrSigmaNoise;
|
---|
[8582] | 84 | public double SigmaNoise {
|
---|
| 85 | get { return Math.Sqrt(sqrSigmaNoise); }
|
---|
| 86 | }
|
---|
[8323] | 87 |
|
---|
| 88 | [Storable]
|
---|
[8982] | 89 | private double[] meanParameter;
|
---|
| 90 | [Storable]
|
---|
| 91 | private double[] covarianceParameter;
|
---|
| 92 |
|
---|
[12819] | 93 | private double[,] l; // used to be storable in previous versions (is calculated lazily now)
|
---|
| 94 | private double[,] x; // scaled training dataset, used to be storable in previous versions (is calculated lazily now)
|
---|
| 95 |
|
---|
| 96 | // BackwardsCompatibility3.4
|
---|
| 97 | #region Backwards compatible code, remove with 3.5
|
---|
| 98 | [Storable(Name = "l")] // restore if available but don't store anymore
|
---|
| 99 | private double[,] l_storable {
|
---|
| 100 | set { this.l = value; }
|
---|
| 101 | get {
|
---|
| 102 | if (trainingDataset == null) return l; // this model has been created with an old version
|
---|
| 103 | else return null; // if the training dataset is available l should not be serialized
|
---|
| 104 | }
|
---|
| 105 | }
|
---|
| 106 | [Storable(Name = "x")] // restore if available but don't store anymore
|
---|
| 107 | private double[,] x_storable {
|
---|
| 108 | set { this.x = value; }
|
---|
| 109 | get {
|
---|
| 110 | if (trainingDataset == null) return x; // this model has been created with an old version
|
---|
| 111 | else return null; // if the training dataset is available x should not be serialized
|
---|
| 112 | }
|
---|
| 113 | }
|
---|
| 114 | #endregion
|
---|
| 115 |
|
---|
| 116 |
|
---|
[8982] | 117 | [Storable]
|
---|
[12819] | 118 | private IDataset trainingDataset; // it is better to store the original training dataset completely because this is more efficient in persistence
|
---|
| 119 | [Storable]
|
---|
| 120 | private int[] trainingRows;
|
---|
[8323] | 121 |
|
---|
| 122 | [Storable]
|
---|
[8463] | 123 | private Scaling inputScaling;
|
---|
[8323] | 124 |
|
---|
| 125 |
|
---|
| 126 | [StorableConstructor]
|
---|
| 127 | private GaussianProcessModel(bool deserializing) : base(deserializing) { }
|
---|
| 128 | private GaussianProcessModel(GaussianProcessModel original, Cloner cloner)
|
---|
| 129 | : base(original, cloner) {
|
---|
| 130 | this.meanFunction = cloner.Clone(original.meanFunction);
|
---|
| 131 | this.covarianceFunction = cloner.Clone(original.covarianceFunction);
|
---|
[13118] | 132 | if (original.inputScaling != null)
|
---|
| 133 | this.inputScaling = cloner.Clone(original.inputScaling);
|
---|
[12819] | 134 | this.trainingDataset = cloner.Clone(original.trainingDataset);
|
---|
[8323] | 135 | this.negativeLogLikelihood = original.negativeLogLikelihood;
|
---|
[14899] | 136 | this.negativeLooPredictiveProbability = original.negativeLooPredictiveProbability;
|
---|
[8416] | 137 | this.sqrSigmaNoise = original.sqrSigmaNoise;
|
---|
[8982] | 138 | if (original.meanParameter != null) {
|
---|
| 139 | this.meanParameter = (double[])original.meanParameter.Clone();
|
---|
| 140 | }
|
---|
| 141 | if (original.covarianceParameter != null) {
|
---|
| 142 | this.covarianceParameter = (double[])original.covarianceParameter.Clone();
|
---|
| 143 | }
|
---|
[8416] | 144 |
|
---|
| 145 | // shallow copies of arrays because they cannot be modified
|
---|
[12819] | 146 | this.trainingRows = original.trainingRows;
|
---|
[8323] | 147 | this.allowedInputVariables = original.allowedInputVariables;
|
---|
| 148 | this.alpha = original.alpha;
|
---|
| 149 | this.l = original.l;
|
---|
| 150 | this.x = original.x;
|
---|
| 151 | }
|
---|
[12509] | 152 | public GaussianProcessModel(IDataset ds, string targetVariable, IEnumerable<string> allowedInputVariables, IEnumerable<int> rows,
|
---|
[13118] | 153 | IEnumerable<double> hyp, IMeanFunction meanFunction, ICovarianceFunction covarianceFunction,
|
---|
| 154 | bool scaleInputs = true)
|
---|
[13941] | 155 | : base(targetVariable) {
|
---|
[8323] | 156 | this.name = ItemName;
|
---|
| 157 | this.description = ItemDescription;
|
---|
[8416] | 158 | this.meanFunction = (IMeanFunction)meanFunction.Clone();
|
---|
| 159 | this.covarianceFunction = (ICovarianceFunction)covarianceFunction.Clone();
|
---|
[8323] | 160 | this.allowedInputVariables = allowedInputVariables.ToArray();
|
---|
| 161 |
|
---|
| 162 |
|
---|
[8416] | 163 | int nVariables = this.allowedInputVariables.Length;
|
---|
[8982] | 164 | meanParameter = hyp
|
---|
[8416] | 165 | .Take(this.meanFunction.GetNumberOfParameters(nVariables))
|
---|
[8982] | 166 | .ToArray();
|
---|
| 167 |
|
---|
| 168 | covarianceParameter = hyp.Skip(this.meanFunction.GetNumberOfParameters(nVariables))
|
---|
| 169 | .Take(this.covarianceFunction.GetNumberOfParameters(nVariables))
|
---|
| 170 | .ToArray();
|
---|
[8473] | 171 | sqrSigmaNoise = Math.Exp(2.0 * hyp.Last());
|
---|
[13160] | 172 | try {
|
---|
| 173 | CalculateModel(ds, rows, scaleInputs);
|
---|
[14843] | 174 | } catch (alglib.alglibexception ae) {
|
---|
[13160] | 175 | // wrap exception so that calling code doesn't have to know about alglib implementation
|
---|
| 176 | throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
|
---|
| 177 | }
|
---|
[8323] | 178 | }
|
---|
| 179 |
|
---|
[13118] | 180 | private void CalculateModel(IDataset ds, IEnumerable<int> rows, bool scaleInputs = true) {
|
---|
[12819] | 181 | this.trainingDataset = (IDataset)ds.Clone();
|
---|
| 182 | this.trainingRows = rows.ToArray();
|
---|
[13118] | 183 | this.inputScaling = scaleInputs ? new Scaling(ds, allowedInputVariables, rows) : null;
|
---|
[8323] | 184 |
|
---|
[13118] | 185 | x = GetData(ds, this.allowedInputVariables, this.trainingRows, this.inputScaling);
|
---|
| 186 |
|
---|
| 187 | IEnumerable<double> y;
|
---|
[13941] | 188 | y = ds.GetDoubleValues(TargetVariable, rows);
|
---|
[13118] | 189 |
|
---|
[8323] | 190 | int n = x.GetLength(0);
|
---|
| 191 |
|
---|
[13721] | 192 | var columns = Enumerable.Range(0, x.GetLength(1)).ToArray();
|
---|
[12819] | 193 | // calculate cholesky decomposed (lower triangular) covariance matrix
|
---|
[13721] | 194 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
|
---|
[12819] | 195 | this.l = CalculateL(x, cov, sqrSigmaNoise);
|
---|
| 196 |
|
---|
| 197 | // calculate mean
|
---|
[13721] | 198 | var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
|
---|
[8982] | 199 | double[] m = Enumerable.Range(0, x.GetLength(0))
|
---|
| 200 | .Select(r => mean.Mean(x, r))
|
---|
| 201 | .ToArray();
|
---|
| 202 |
|
---|
[8323] | 203 | // calculate sum of diagonal elements for likelihood
|
---|
| 204 | double diagSum = Enumerable.Range(0, n).Select(i => Math.Log(l[i, i])).Sum();
|
---|
| 205 |
|
---|
| 206 | // solve for alpha
|
---|
| 207 | double[] ym = y.Zip(m, (a, b) => a - b).ToArray();
|
---|
| 208 |
|
---|
[12819] | 209 | int info;
|
---|
| 210 | alglib.densesolverreport denseSolveRep;
|
---|
| 211 |
|
---|
[8323] | 212 | alglib.spdmatrixcholeskysolve(l, n, false, ym, out info, out denseSolveRep, out alpha);
|
---|
| 213 | for (int i = 0; i < alpha.Length; i++)
|
---|
| 214 | alpha[i] = alpha[i] / sqrSigmaNoise;
|
---|
| 215 | negativeLogLikelihood = 0.5 * Util.ScalarProd(ym, alpha) + diagSum + (n / 2.0) * Math.Log(2.0 * Math.PI * sqrSigmaNoise);
|
---|
| 216 |
|
---|
| 217 | // derivatives
|
---|
| 218 | int nAllowedVariables = x.GetLength(1);
|
---|
| 219 |
|
---|
[8463] | 220 | alglib.matinvreport matInvRep;
|
---|
[8475] | 221 | double[,] lCopy = new double[l.GetLength(0), l.GetLength(1)];
|
---|
| 222 | Array.Copy(l, lCopy, lCopy.Length);
|
---|
[8323] | 223 |
|
---|
[8475] | 224 | alglib.spdmatrixcholeskyinverse(ref lCopy, n, false, out info, out matInvRep);
|
---|
[8463] | 225 | if (info != 1) throw new ArgumentException("Can't invert matrix to calculate gradients.");
|
---|
[14899] | 226 |
|
---|
| 227 | // LOOCV log predictive probability (GPML page 116 and 117)
|
---|
| 228 | var sumLoo = 0.0;
|
---|
| 229 | var ki = new double[n];
|
---|
[8323] | 230 | for (int i = 0; i < n; i++) {
|
---|
[14899] | 231 | for (int j = 0; j < n; j++) ki[j] = cov.Covariance(x, i, j);
|
---|
| 232 | var yi = Util.ScalarProd(ki, alpha);
|
---|
| 233 | var yi_loo = yi - alpha[i] / lCopy[i, i] / sqrSigmaNoise;
|
---|
| 234 | var s2_loo = sqrSigmaNoise / lCopy[i, i];
|
---|
| 235 | var err = ym[i] - yi_loo;
|
---|
| 236 | var nll_loo = Math.Log(s2_loo) + err * err / s2_loo;
|
---|
| 237 | sumLoo += nll_loo;
|
---|
| 238 | }
|
---|
| 239 | sumLoo += n * Math.Log(2 * Math.PI);
|
---|
| 240 | negativeLooPredictiveProbability = 0.5 * sumLoo;
|
---|
| 241 |
|
---|
| 242 | for (int i = 0; i < n; i++) {
|
---|
[8463] | 243 | for (int j = 0; j <= i; j++)
|
---|
[8475] | 244 | lCopy[i, j] = lCopy[i, j] / sqrSigmaNoise - alpha[i] * alpha[j];
|
---|
[8323] | 245 | }
|
---|
| 246 |
|
---|
[8475] | 247 | double noiseGradient = sqrSigmaNoise * Enumerable.Range(0, n).Select(i => lCopy[i, i]).Sum();
|
---|
[8323] | 248 |
|
---|
| 249 | double[] meanGradients = new double[meanFunction.GetNumberOfParameters(nAllowedVariables)];
|
---|
[8982] | 250 | for (int k = 0; k < meanGradients.Length; k++) {
|
---|
[13721] | 251 | var meanGrad = new double[alpha.Length];
|
---|
| 252 | for (int g = 0; g < meanGrad.Length; g++)
|
---|
| 253 | meanGrad[g] = mean.Gradient(x, g, k);
|
---|
[8982] | 254 | meanGradients[k] = -Util.ScalarProd(meanGrad, alpha);
|
---|
[8323] | 255 | }
|
---|
| 256 |
|
---|
| 257 | double[] covGradients = new double[covarianceFunction.GetNumberOfParameters(nAllowedVariables)];
|
---|
[8366] | 258 | if (covGradients.Length > 0) {
|
---|
| 259 | for (int i = 0; i < n; i++) {
|
---|
[8484] | 260 | for (int j = 0; j < i; j++) {
|
---|
[13784] | 261 | var g = cov.CovarianceGradient(x, i, j);
|
---|
[8484] | 262 | for (int k = 0; k < covGradients.Length; k++) {
|
---|
| 263 | covGradients[k] += lCopy[i, j] * g[k];
|
---|
[8366] | 264 | }
|
---|
[8323] | 265 | }
|
---|
[8484] | 266 |
|
---|
[13784] | 267 | var gDiag = cov.CovarianceGradient(x, i, i);
|
---|
[8484] | 268 | for (int k = 0; k < covGradients.Length; k++) {
|
---|
| 269 | // diag
|
---|
| 270 | covGradients[k] += 0.5 * lCopy[i, i] * gDiag[k];
|
---|
| 271 | }
|
---|
[8323] | 272 | }
|
---|
| 273 | }
|
---|
| 274 |
|
---|
[8484] | 275 | hyperparameterGradients =
|
---|
[8473] | 276 | meanGradients
|
---|
| 277 | .Concat(covGradients)
|
---|
| 278 | .Concat(new double[] { noiseGradient }).ToArray();
|
---|
[8484] | 279 |
|
---|
[8323] | 280 | }
|
---|
| 281 |
|
---|
[13118] | 282 | private static double[,] GetData(IDataset ds, IEnumerable<string> allowedInputs, IEnumerable<int> rows, Scaling scaling) {
|
---|
| 283 | if (scaling != null) {
|
---|
[14854] | 284 | // BackwardsCompatibility3.3
|
---|
| 285 | #region Backwards compatible code, remove with 3.4
|
---|
[14843] | 286 | // TODO: completely remove Scaling class
|
---|
[14854] | 287 | List<string> variablesList = allowedInputs.ToList();
|
---|
| 288 | List<int> rowsList = rows.ToList();
|
---|
[14843] | 289 |
|
---|
[14854] | 290 | double[,] matrix = new double[rowsList.Count, variablesList.Count];
|
---|
| 291 |
|
---|
| 292 | int col = 0;
|
---|
| 293 | foreach (string column in variablesList) {
|
---|
| 294 | var values = scaling.GetScaledValues(ds, column, rowsList);
|
---|
| 295 | int row = 0;
|
---|
| 296 | foreach (var value in values) {
|
---|
| 297 | matrix[row, col] = value;
|
---|
| 298 | row++;
|
---|
| 299 | }
|
---|
| 300 | col++;
|
---|
[14843] | 301 | }
|
---|
[14854] | 302 | return matrix;
|
---|
| 303 | #endregion
|
---|
[13118] | 304 | } else {
|
---|
[14843] | 305 | return ds.ToArray(allowedInputs, rows);
|
---|
[13118] | 306 | }
|
---|
[12819] | 307 | }
|
---|
[8323] | 308 |
|
---|
[12819] | 309 | private static double[,] CalculateL(double[,] x, ParameterizedCovarianceFunction cov, double sqrSigmaNoise) {
|
---|
| 310 | int n = x.GetLength(0);
|
---|
| 311 | var l = new double[n, n];
|
---|
| 312 |
|
---|
| 313 | // calculate covariances
|
---|
| 314 | for (int i = 0; i < n; i++) {
|
---|
| 315 | for (int j = i; j < n; j++) {
|
---|
| 316 | l[j, i] = cov.Covariance(x, i, j) / sqrSigmaNoise;
|
---|
| 317 | if (j == i) l[j, i] += 1.0;
|
---|
| 318 | }
|
---|
| 319 | }
|
---|
| 320 |
|
---|
| 321 | // cholesky decomposition
|
---|
| 322 | var res = alglib.trfac.spdmatrixcholesky(ref l, n, false);
|
---|
| 323 | if (!res) throw new ArgumentException("Matrix is not positive semidefinite");
|
---|
| 324 | return l;
|
---|
| 325 | }
|
---|
| 326 |
|
---|
| 327 |
|
---|
[8323] | 328 | public override IDeepCloneable Clone(Cloner cloner) {
|
---|
| 329 | return new GaussianProcessModel(this, cloner);
|
---|
| 330 | }
|
---|
| 331 |
|
---|
[8982] | 332 | // is called by the solution creator to set all parameter values of the covariance and mean function
|
---|
| 333 | // to the optimized values (necessary to make the values visible in the GUI)
|
---|
| 334 | public void FixParameters() {
|
---|
| 335 | covarianceFunction.SetParameter(covarianceParameter);
|
---|
| 336 | meanFunction.SetParameter(meanParameter);
|
---|
| 337 | covarianceParameter = new double[0];
|
---|
| 338 | meanParameter = new double[0];
|
---|
| 339 | }
|
---|
| 340 |
|
---|
[8323] | 341 | #region IRegressionModel Members
|
---|
[13941] | 342 | public override IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
|
---|
[8323] | 343 | return GetEstimatedValuesHelper(dataset, rows);
|
---|
| 344 | }
|
---|
[13941] | 345 | public override IRegressionSolution CreateRegressionSolution(IRegressionProblemData problemData) {
|
---|
[8528] | 346 | return new GaussianProcessRegressionSolution(this, new RegressionProblemData(problemData));
|
---|
[8323] | 347 | }
|
---|
| 348 | #endregion
|
---|
| 349 |
|
---|
[8623] | 350 |
|
---|
[12509] | 351 | private IEnumerable<double> GetEstimatedValuesHelper(IDataset dataset, IEnumerable<int> rows) {
|
---|
[13160] | 352 | try {
|
---|
| 353 | if (x == null) {
|
---|
| 354 | x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
|
---|
| 355 | }
|
---|
| 356 | int n = x.GetLength(0);
|
---|
[12819] | 357 |
|
---|
[13160] | 358 | double[,] newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
|
---|
| 359 | int newN = newX.GetLength(0);
|
---|
[12819] | 360 |
|
---|
[13721] | 361 | var Ks = new double[newN][];
|
---|
| 362 | var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
|
---|
| 363 | var mean = meanFunction.GetParameterizedMeanFunction(meanParameter, columns);
|
---|
[13160] | 364 | var ms = Enumerable.Range(0, newX.GetLength(0))
|
---|
| 365 | .Select(r => mean.Mean(newX, r))
|
---|
| 366 | .ToArray();
|
---|
[13721] | 367 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
|
---|
[13160] | 368 | for (int i = 0; i < newN; i++) {
|
---|
[13721] | 369 | Ks[i] = new double[n];
|
---|
[13160] | 370 | for (int j = 0; j < n; j++) {
|
---|
[13721] | 371 | Ks[i][j] = cov.CrossCovariance(x, newX, j, i);
|
---|
[13160] | 372 | }
|
---|
[8323] | 373 | }
|
---|
[13160] | 374 |
|
---|
| 375 | return Enumerable.Range(0, newN)
|
---|
[13721] | 376 | .Select(i => ms[i] + Util.ScalarProd(Ks[i], alpha));
|
---|
[14843] | 377 | } catch (alglib.alglibexception ae) {
|
---|
[13160] | 378 | // wrap exception so that calling code doesn't have to know about alglib implementation
|
---|
| 379 | throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
|
---|
[8323] | 380 | }
|
---|
| 381 | }
|
---|
[8473] | 382 |
|
---|
[14095] | 383 | public IEnumerable<double> GetEstimatedVariances(IDataset dataset, IEnumerable<int> rows) {
|
---|
[13160] | 384 | try {
|
---|
| 385 | if (x == null) {
|
---|
| 386 | x = GetData(trainingDataset, allowedInputVariables, trainingRows, inputScaling);
|
---|
| 387 | }
|
---|
| 388 | int n = x.GetLength(0);
|
---|
[12819] | 389 |
|
---|
[13160] | 390 | var newX = GetData(dataset, allowedInputVariables, rows, inputScaling);
|
---|
| 391 | int newN = newX.GetLength(0);
|
---|
[8473] | 392 |
|
---|
[13160] | 393 | var kss = new double[newN];
|
---|
| 394 | double[,] sWKs = new double[n, newN];
|
---|
[13721] | 395 | var columns = Enumerable.Range(0, newX.GetLength(1)).ToArray();
|
---|
| 396 | var cov = covarianceFunction.GetParameterizedCovarianceFunction(covarianceParameter, columns);
|
---|
[8473] | 397 |
|
---|
[13160] | 398 | if (l == null) {
|
---|
| 399 | l = CalculateL(x, cov, sqrSigmaNoise);
|
---|
| 400 | }
|
---|
[12819] | 401 |
|
---|
[13160] | 402 | // for stddev
|
---|
| 403 | for (int i = 0; i < newN; i++)
|
---|
| 404 | kss[i] = cov.Covariance(newX, i, i);
|
---|
[8473] | 405 |
|
---|
[13160] | 406 | for (int i = 0; i < newN; i++) {
|
---|
| 407 | for (int j = 0; j < n; j++) {
|
---|
| 408 | sWKs[j, i] = cov.CrossCovariance(x, newX, j, i) / Math.Sqrt(sqrSigmaNoise);
|
---|
| 409 | }
|
---|
[8473] | 410 | }
|
---|
| 411 |
|
---|
[13160] | 412 | // for stddev
|
---|
| 413 | alglib.ablas.rmatrixlefttrsm(n, newN, l, 0, 0, false, false, 0, ref sWKs, 0, 0);
|
---|
[8473] | 414 |
|
---|
[13160] | 415 | for (int i = 0; i < newN; i++) {
|
---|
[13721] | 416 | var col = Util.GetCol(sWKs, i).ToArray();
|
---|
| 417 | var sumV = Util.ScalarProd(col, col);
|
---|
[13160] | 418 | kss[i] += sqrSigmaNoise; // kss is V(f), add noise variance of predictive distibution to get V(y)
|
---|
| 419 | kss[i] -= sumV;
|
---|
| 420 | if (kss[i] < 0) kss[i] = 0;
|
---|
| 421 | }
|
---|
| 422 | return kss;
|
---|
[14843] | 423 | } catch (alglib.alglibexception ae) {
|
---|
[13160] | 424 | // wrap exception so that calling code doesn't have to know about alglib implementation
|
---|
| 425 | throw new ArgumentException("There was a problem in the calculation of the Gaussian process model", ae);
|
---|
[8473] | 426 | }
|
---|
| 427 | }
|
---|
[13921] | 428 |
|
---|
[8323] | 429 | }
|
---|
| 430 | }
|
---|