Commits

Ruben Martinez-Cantin committed 28b6f89

Full Bayesian approach fully working. Solved bug in kernel parameters prior.

Comments (0)

Files changed (6)

 Start by reading the \ref install and the \ref reference
 
 You can also find more details at:
-<http://bitbucket.org/rmcantin/bayesian-optimization/wiki/Home>
+<http://bitbucket.org/rmcantin/bayesopt/wiki/Home>
 
 **Important:** This code is free to use. However, if you are using,
 or plan to use, the library, specially if it is for research or

app/bo_display.cpp

 {
   size_t dim = 1;
   bopt_params parameters = initialize_parameters_to_default();
-  parameters.n_init_samples = 7;
+  parameters.n_init_samples = 10;
   parameters.n_iter_relearn = 0;
   parameters.n_iterations = 150;
   parameters.surr_name = S_STUDENT_T_PROCESS_NORMAL_INV_GAMMA;
   parameters.kernel.hp_mean[0] = 1;
-  parameters.kernel.hp_std[0] = 0.1;
+  parameters.kernel.hp_std[0] = 5;
   parameters.kernel.n_hp = 1;
+  //  parameters.noise = 0.01;
+  // parameters.l_type = L_ML;
   // parameters.mean.name = "mZero";
   //  parameters.crit_name = "cHedge(cEI,cLCB,cExpReturn,cOptimisticSampling)";
   // parameters.epsilon = 0.0;

include/parameters.h

 
   /* Nonparametric process "parameters" */
   const double KERNEL_THETA    = 1.0;
-  const double KERNEL_SIGMA    = 10.0;
+  const double KERNEL_SIGMA    = 100.0;
   const double MEAN_MU         = 1.0;
   const double MEAN_SIGMA      = 1000.0;
   const double PRIOR_ALPHA     = 1.0;

src/inneroptimization.cpp

 
     nlopt_opt opt;
     double (*fpointer)(unsigned int, const double *, double *, void *);
-    double coef = 1.0;  //Percentaje of resources used in local optimization
+    double coef = 0.8;  //Percentaje of resources used in local optimization
 
     /* algorithm and dims */
     if (alg == LBFGS)                                     //Require gradient

src/nonparametricprocess.cpp

 
 #include "gaussian_process.hpp"
 #include "gaussian_process_ml.hpp"
-//#include "gaussian_process_ign.hpp"
 #include "student_t_process_jef.hpp"
 #include "student_t_process_nig.hpp"
 
     mSigma(parameters.sigma_s), dim_(dim)
   { 
     mMinIndex = 0;     mMaxIndex = 0;   
-    setAlgorithm(BOBYQA);
-    setLimits(0.,100.);
+    if (parameters.l_type == L_ML)
+      {
+	setAlgorithm(BOBYQA);    // local search to avoid underfitting
+      }
+    else
+      {
+	setAlgorithm(COMBINED);
+      }
+    setLimits(1e-10,100.);
     setLearnType(parameters.l_type);
     setKernel(parameters.kernel,dim);
     setMean(parameters.mean,dim);
       {
 	if (priorKernel[i].standard_deviation() > 0)
 	  {
-	    prior += log(boost::math::pdf(priorKernel[i],th(i)));
+	    prior -= log(boost::math::pdf(priorKernel[i],th(i)));
 	  }
       }
     return prior;

src/student_t_process_nig.cpp

 ------------------------------------------------------------------------
 */
 
+#include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/numeric/ublas/banded.hpp>
 #include "log.hpp"
 #include "student_t_process_nig.hpp"
     double sPred = sqrt( mSigmaMap * (kq - inner_prod(v,v) 
 				   + inner_prod(rho,rho)));
 
+    if ((boost::math::isnan(yPred)) || (boost::math::isnan(sPred)))
+      {
+	FILE_LOG(logERROR) << "Error in prediction. NaN found.";
+	throw 1;
+      }
+					
+
     d_->setMeanAndStd(yPred,sPred);
     return d_;
   }
 
   double StudentTProcessNIG::negativeLogLikelihood()
   {
-    /* TODO: For testing */
-    return negativeTotalLogLikelihood();
+    matrixd KK = computeCorrMatrix();
+    const size_t n = KK.size1();
+    const size_t p = mMean->nFeatures();
+    const size_t nalpha = (n+2*mAlpha);
+
+    vectord v0 = mGPY - prod(trans(mFeatM),mW0);
+    matrixd WW = zmatrixd(p,p);  //TODO: diagonal matrix
+    utils::addToDiagonal(WW,mInvVarW);
+    matrixd FW = prod(trans(mFeatM),WW);
+    KK += prod(FW,mFeatM);
+    matrixd BB(n,n);
+    utils::cholesky_decompose(KK,BB);
+    inplace_solve(BB,v0,ublas::lower_tag());
+    double zz = inner_prod(v0,v0);
+    double sigmaMap = (mBeta/mAlpha + zz)/nalpha;
+
+    double lik = nalpha/2 * log(1+zz/(2*mBeta*sigmaMap));
+    lik += utils::log_trace(BB);
+    lik += n/2 * log(sigmaMap);
+    return lik;
   }
 
 
 
     vectord v0 = mGPY - prod(trans(mFeatM),mW0);
     //TODO: check for "cheaper" version
-    matrixd KK = prod(mL,trans(mL));
+    //matrixd KK = prod(mL,trans(mL));
+    matrixd KK = computeCorrMatrix();
     matrixd WW = zmatrixd(p,p);  //TODO: diagonal matrix
     utils::addToDiagonal(WW,mInvVarW);
     matrixd FW = prod(trans(mFeatM),WW);
     
     int dof = static_cast<int>(n+2*mAlpha);
     
+    if ((boost::math::isnan(mWMap(0))) || (boost::math::isnan(mSigmaMap)))
+      {
+	FILE_LOG(logERROR) << "Error in precomputed prediction. NaN found.";
+	throw 1;
+      }
+
+
     if (dof <= 0)  
       {
 	FILE_LOG(logERROR) << "ERROR: Incorrect alpha. Dof invalid.";