I think your package could be very handy for my problem; thanks for developing it.
I prototyped some R code, which works fine:
library(bigsnpr)
bigsnp <- snp_attachExtdata()
G <- bigsnp$genotypes
simu <- snp_simuPheno(G, 0.2, 500, alpha = 0.5)
log_var <- log(big_colstats(G, ind.col = simu$set)$var)
beta2 <- simu$effects^2
FUN <- function(x, log_var, beta2) {
S <- 1 + x[[1]]; sigma2 <- x[[2]]
S * sum(log_var) + length(log_var) * log(sigma2) + sum(beta2 / exp(S * log_var)) / sigma2
}
DER <- function(x, log_var, beta2) {
S <- 1 + x[[1]]; sigma2 <- x[[2]]
res1 <- sum(log_var) - sum(log_var * beta2 / exp(S * log_var)) / sigma2
res2 <- length(log_var) / sigma2 - sum(beta2 / exp(S * log_var)) / sigma2^2
c(res1, res2)
}
optim(par = c(-0.5, 0.2 / 500), fn = FUN, method = "L-BFGS-B", gr = DER,
lower = c(-2, 0.2 / 5000), upper = c(1, 0.2 / 50),
log_var = log_var, beta2 = beta2)
Then I tried implementing it in C++ using your package:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
#include <roptim.h>
// [[Rcpp::depends(roptim)]]
using namespace roptim;
class MLE : public Functor {
public:
MLE(const arma::vec& a_, const arma::vec& b_) : a(a_), b(b_) {}
double operator()(const arma::vec& par) override {
double S = par[0] + 1;
double sigma2 = par[1];
int m = a.size();
double sum_a = 0, sum_c = 0;
for (int j = 0; j < m; j++) {
sum_a += a[j];
sum_c += b[j] / ::exp(S * a[j]);
}
return S * sum_a + m * ::log(sigma2) + sum_c / sigma2;
}
void Gradient(const arma::vec& par, arma::vec& gr) override {
double S = par[0] + 1;
double sigma2 = par[1];
int m = a.size();
double sum_a = 0, sum_c = 0, sum_ac = 0;
for (int j = 0; j < m; j++) {
double c_j = b[j] / ::exp(S * a[j]);
sum_a += a[j];
sum_c += c_j;
sum_ac += a[j] * c_j;
}
gr[0] = sum_a - sum_ac / sigma2;
gr[1] = (m - sum_c / sigma2) / sigma2;
}
private:
arma::vec a;
arma::vec b;
};
// [[Rcpp::export]]
arma::vec test_MLE(const arma::vec& a_, const arma::vec& b_,
arma::vec par, const arma::vec& lower, const arma::vec& upper) {
MLE mle(a_, b_);
Roptim<MLE> opt("L-BFGS-B");
opt.set_lower(lower);
opt.set_upper(upper);
opt.set_hessian(false);
opt.control.maxit = 100;
opt.minimize(mle, par);
Rcpp::Rcout << "-------------------------" << std::endl;
opt.print();
return par;
}
Then in R:
Rcpp::sourceCpp('tmp-tests/proto-MLE-optim-C++.cpp')
x <- c(0.5, 0.2 / 500)
test_MLE(log_var, beta2, x,
lower = c(-2, 0.2 / 5000), upper = c(1, 0.2 / 50))
but it crashes..
I have been trying to debug it, but I am not sure if it comes from a silly mistake on my part, or something I didn't understand about how to use your package.
Any idea?