estimate_alpha.h

Go to the documentation of this file.

 
#pragma once
 
#include "ols.h"
 
namespace slope {
 
template<typename T>
double
estimateNoise(Eigen::EigenBase<T>& 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] =
      detail::fitOls(x.derived(), y, fit_intercept);
    residuals = y - x.derived() * 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);
}
 
} // namespace slope