24 #include <boost/math/special_functions/fpclassify.hpp> 25 #include <boost/numeric/ublas/banded.hpp> 26 #include "ublas_trace.hpp" 34 namespace ublas = boost::numeric::ublas;
36 GaussianProcessNormal::GaussianProcessNormal(
size_t dim,
41 HierarchicalGaussianProcess(dim,params,data,mean,eng),
42 mW0(params.mean.coef_mean.size()), mInvVarW(params.mean.coef_mean.size()),
43 mD(params.mean.coef_mean.size(),params.mean.coef_mean.size())
46 mW0 = params.mean.coef_mean;
47 for (
size_t ii = 0; ii < params.mean.coef_std.size(); ++ii)
49 double varii = params.mean.coef_std[ii] * params.mean.coef_std[ii];
50 mInvVarW(ii) = 1/varii;
52 d_ =
new GaussianDistribution(eng);
57 GaussianProcessNormal::~GaussianProcessNormal()
63 ProbabilityDistribution*
66 const double kq = computeSelfCorrelation(query);
67 const vectord phi = mMean.getFeatures(query);
69 vectord v = computeCrossCorrelation(query);
71 inplace_solve(
mL,v,ublas::lower_tag());
73 vectord rq = phi - prod(v,mKF);
76 inplace_solve(mD,rho,ublas::lower_tag());
78 double yPred = inner_prod(phi,mWMap) + inner_prod(v,mVf);
79 double sPred = sqrt(
mSigma * (kq - inner_prod(v,v)
80 + inner_prod(rho,rho)));
82 if ((boost::math::isnan(yPred)) || (boost::math::isnan(sPred)))
84 throw std::runtime_error(
"Error in prediction. NaN found.");
96 const size_t n = KK.size1();
97 const size_t p = mMean.nFeatures();
99 vectord v0 = mData.
mY - prod(trans(mMean.
mFeatM),mW0);
100 matrixd WW = zmatrixd(p,p);
101 utils::add_to_diagonal(WW,mInvVarW);
102 matrixd FW = prod(trans(mMean.
mFeatM),WW);
103 KK += prod(FW,mMean.
mFeatM);
106 inplace_solve(BB,v0,ublas::lower_tag());
107 double zz = inner_prod(v0,v0);
109 double lik = 1/(2*
mSigma) * zz;
110 lik += utils::log_trace(BB);
118 const size_t p = mMean.nFeatures();
120 mKF = trans(mMean.
mFeatM);
121 inplace_solve(
mL,mKF,ublas::lower_tag());
124 DD = prod(trans(mKF),mKF);
125 utils::add_to_diagonal(DD,mInvVarW);
128 vectord vn = mData.
mY;
129 inplace_solve(
mL,vn,ublas::lower_tag());
133 mVf = mData.
mY - prod(trans(mMean.
mFeatM),mWMap);
134 inplace_solve(
mL,mVf,ublas::lower_tag());
136 if (boost::math::isnan(mWMap(0)))
138 throw std::runtime_error(
"Error in precomputed prediction. NaN found.");
ProbabilityDistribution * prediction(const vectord &query)
Function that returns the prediction of the GP for a query point in the hypercube [0...
matrixd computeCorrMatrix()
Computes the Correlation (Kernel or Gram) matrix.
void setMeanAndStd(double mean, double std)
Sets the mean and std of the distribution.
void precomputePrediction()
Precompute some values of the prediction that do not depends on the query.
Namespace of the library interface.
Elementwise operations for ublas vector/matrix.
matrixd mL
Cholesky decomposition of the Correlation matrix.
GaussianDistribution * d_
Predictive distributions.
Gaussian probability distribution.
v1 ublas_elementwise_prod(const v1 &a, const v2 &b)
Computes the elementwise product of two vectors or matrices.
double negativeLogLikelihood()
Computes the negative log likelihood and its gradient of the data.
double mSigma
Signal variance.
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
void cholesky_solve(const TRIA &L, MATRIX &x, ublas::lower)
solve system L L^T x = b inplace
matrixd mFeatM
Value of the mean features at the input points.
Gaussian process with normal prior on the parameters.