estimate_alpha.h

Go to the documentation of this file.

 
#pragma once
 
#include "logger.h"
#include "ols.h"
#include "slope.h"
#include "utils.h"
 
namespace slope {
 
template<typename MatrixType>
double
estimateNoise(MatrixType& x, Eigen::MatrixXd& y, const bool fit_intercept)
{
  int n = x.rows();
  int p = x.cols();
 
  Eigen::VectorXd residuals;
  int df; // Degrees of freedom
 
  if (p == 0) {
    // Handle the case when X has zero columns
    residuals = y;
    df = n - 1;
 
    if (fit_intercept) {
      y.array() -= y.mean();
      df -= 1;
    }
  } else {
    // Normal case with predictors
    auto [ols_intercept, ols_coefs] = fitOls(x, y, fit_intercept);
    residuals = y - x * ols_coefs;
 
    if (fit_intercept) {
      residuals.array() -= ols_intercept;
    }
 
    df = n - p - static_cast<int>(fit_intercept) - 1;
  }
 
  assert(df > 0);
 
  return residuals.norm() / std::sqrt(df);
}
 
template<typename MatrixType>
SlopePath
estimateAlpha(MatrixType& x, Eigen::MatrixXd& y, const Slope& model)
{
  int n = x.rows();
  int p = x.cols();
 
  int alpha_est_maxit = model.getAlphaEstimationMaxIterations();
 
  // TODO: Possibly we could just run the path instead. We could use the
  // selected predictors at its end as the selected set and run these in an OLS.
  // There would no longer be any problem with cyclic behavior of this
  // algorithm, but on the other hand the model is likely very unstable.
  std::vector<int> selected;
 
  Eigen::ArrayXd alpha(1);
 
  Slope model_copy = model;
 
  model_copy.setAlphaType("path"); // Otherwise we would be in an infinite loop
 
  SlopePath result;
 
  // Estimate the noise level, if possible
  if (n >= p + 30) {
    alpha(0) = estimateNoise(x, y, model_copy.getFitIntercept()) / n;
    result = model_copy.path(x, y, alpha);
  } else {
    for (int it = 0; it < alpha_est_maxit; ++it) {
 
      auto x_selected = subsetCols(x, selected);
 
      std::vector<int> selected_prev = selected;
      selected.clear();
 
      alpha(0) = estimateNoise(x_selected, y, model_copy.getFitIntercept()) / n;
 
      // TODO: If changes in alpha are small between two steps, then it should
      // be easy to screen, but we are not using that possibility here.
      result = model_copy.path(x, y, alpha);
      auto coefs = result.getCoefs().back();
 
      for (typename Eigen::SparseMatrix<double>::InnerIterator it(coefs, 0); it;
           ++it) {
        selected.emplace_back(it.row());
      }
 
      if (selected == selected_prev) {
        return result;
      }
 
      if (static_cast<int>(selected.size()) >=
          n + model_copy.getFitIntercept()) {
        throw std::runtime_error(
          "selected >= n - 1 variables, cannot estimate variance");
      }
    }
 
    slope::WarningLogger::addWarning(
      slope::WarningCode::MAXIT_REACHED,
      "Maximum iterations reached in alpha estimation");
  }
 
  return result;
}
} // namespace slope