# Bayesian-Optimization / src / basicgaussprocess.cpp

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117``` ```#include #include // for normal_distribution #include "basicgaussprocess.hpp" #include "cholesky.hpp" #include "trace.hpp" BasicGaussianProcess::BasicGaussianProcess( double theta, double noise): NonParametricProcess(theta,noise) { setAlgorithm(bobyqa); setLimits(0.,100.); } // Constructor BasicGaussianProcess::~BasicGaussianProcess() {} // Default destructor double BasicGaussianProcess::negativeLogLikelihood(double &grad, size_t index) { matrixd K = computeCorrMatrix(mRegularizer,0); size_t n = K.size1(); matrixd L(n,n); cholesky_decompose(K,L); // Compute the likelihood vectord alpha(mGPY); boost::numeric::ublas::inplace_solve(L,alpha,boost::numeric::ublas::lower_tag()); double loglik = .5*inner_prod(mGPY,alpha) + trace(L) + n*0.91893853320467; //log(2*pi)/2 #if 0 // BUG: // Compute the ith derivative if (index > 0) { matrixd inverse = eyed(n); boost::numeric::ublas::inplace_solve(L,inverse,lower_tag()); matrixd dK = computeCorrMatrix(0,index); grad = -.5 * trace_prod(outer_prod(alpha,alpha) - inverse, dK); } #endif return loglik; } int BasicGaussianProcess::prediction( const vectord &query, double& yPred, double& sPred) { vectord rInvR(mGPXX.size()); double kn; double rInvRr; vectord colR = computeCrossCorrelation(query); kn = correlationFunction(query, query); noalias(rInvR) = prod(colR,mInvR); rInvRr = inner_prod(rInvR,colR); yPred = inner_prod( rInvR, mGPY ); sPred = sqrt(kn - rInvRr); return 1; } double BasicGaussianProcess::negativeExpectedImprovement(double yPred, double sPred, double yMin, size_t g) { using boost::math::factorial; boost::math::normal d; double yDiff = yMin - yPred; double yNorm = yDiff / sPred; if (g == 1) return -1.0 * ( yDiff * cdf(d,yNorm) + sPred * pdf(d,yNorm) ); else { double pdfD = pdf(d,yNorm); double Tm2 = cdf(d,yNorm); double Tm1 = pdfD; double fg = factorial(g); double Tact; double sumEI = pow(yNorm,g)*Tm2 - g*pow(yNorm,g-1)*Tm1; for (unsigned int ii = 2; ii < g; ii++) { Tact = (ii-1)*Tm2 - pdfD*pow(yNorm,ii-1); sumEI += pow(-1.0,ii)* (fg / ( factorial(ii)*factorial(g-ii) ) )* pow(yNorm,g-ii)*Tact; //roll-up Tm2 = Tm1; Tm1 = Tact; } return -1.0 * pow(sPred,g) * sumEI; } } // negativeExpectedImprovement double BasicGaussianProcess::lowerConfidenceBound(double yPred, double sPred, double beta) { return yPred - beta*sPred;; } // lowerConfidenceBound double BasicGaussianProcess::negativeProbabilityOfImprovement(double yPred, double sPred, double yMin, double epsilon) { boost::math::normal d; return -cdf(d,(yMin - yPred + epsilon)/sPred); } // negativeProbabilityOfImprovement ```
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.