[EDIT]
Я ищу способ использования функции optim_lbfgs в Rcppnumeric и RcppEigen с Rcpparmadillo.Я следую пути в Интеграция с Rcpp для библиотек числовых вычислений , но он не работал с ошибкой, cannot declare variable 'obj' to be of abstract type 'scoreftn_mns'.Но теперь я исправил некоторые коды, чтобы он работал.Я определяю beta как Eigen::VectorXd beta(p); как в Rcppnumeric и преобразовываю его в arma :: vec в шаблоне (?).
Вот мой код, который я пытаюсь сделать,
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <RcppNumerical.h>
using namespace Numer;
using namespace arma;
using namespace Rcpp;
// Eigen::Ref<Eigen::VectorXd>
// Eigen::Ref<const Eigen::VectorXd>
class socreftn_mns: public MFuncGrad
{
private:
const arma::vec TIME;
const arma::vec DELTA;
const arma::mat COVARI;
const arma::vec TARGETVEC;
public:
socreftn_mns(const arma::vec Time, const arma::vec Delta, const arma::mat Covari,
const arma::vec targetvec) : TIME(Time), DELTA(Delta), COVARI(Covari), TARGETVEC(targetvec) {}
double f_grad(Constvec& beta, Refvec grad){
arma::vec b_s = arma::vec(beta.data(),beta.size());
int n = COVARI.n_rows;
int p = COVARI.n_cols;
arma::vec zero_vec_p = zeros(p);
arma::mat zero_mat_np = zeros(n,p);
arma::vec tempvec_p(p);
arma::mat tempmat_np(n,p);
arma::vec resid = log(TIME) + COVARI*b_s;
arma::uvec index_resid = sort_index(resid);
TIME(index_resid);
DELTA(index_resid);
COVARI.rows(index_resid);
resid(index_resid);
tempmat_np = zero_mat_np; arma::vec U_inf = zero_vec_p;
for(int it=0; it<n; it++){
tempmat_np = COVARI.row(it) - COVARI.each_row();
U_inf += sum(tempmat_np.each_col()%conv_to<vec>::from((resid>=resid(it))),0).t()*DELTA(it);
}
U_inf = U_inf/n - TARGETVEC;
double objvalue = conv_to<double>::from(sum(pow(U_inf,2)));
double h = 1e-4;
for(int itt=0; itt<p; itt++){
tempvec_p = b_s;
tempvec_p(itt) = tempvec_p(itt) + h;
arma::vec resid_g = log(TIME) + COVARI*tempvec_p;
arma::uvec index_resid_g = sort_index(resid_g);
TIME(index_resid_g);
DELTA(index_resid_g);
COVARI.rows(index_resid_g);
resid(index_resid_g);
tempmat_np = zero_mat_np; arma::vec score_grad = zero_vec_p;
for(int it=0; it<n; it++){
tempmat_np = COVARI.row(it) - COVARI.each_row();
score_grad += sum(tempmat_np.each_col()%conv_to<vec>::from((resid_g>=resid_g(it))),0).t()*DELTA(it);
}
score_grad = score_grad/n - TARGETVEC;
double score_objvalue = conv_to<double>::from(sum(pow(score_grad,2)));
grad(itt) = (score_objvalue-objvalue)/h;
}
// beta = Eigen::Ref(b_s);
return objvalue;
}
};
// [[Rcpp::export]]
Rcpp::NumericVector aftsrr_bfgs(arma::vec Time, arma::vec Delta, arma::mat Covari, arma::vec targetvec){
const arma::vec TIME = Time;
const arma::vec DELTA = Delta;
const arma::mat COVARI = Covari;
const arma::vec TARGETVEC = targetvec;
int p = COVARI.n_cols;
// Score Function
socreftn_mns obj(TIME, DELTA, COVARI, TARGETVEC);
// Initial Guess
Eigen::VectorXd beta(p);
beta.setOnes();
double fopt;
int status = optim_lbfgs(obj, beta, fopt);
if(status < 0)
Rcpp::stop("fail to converge");
return Rcpp::wrap(beta);
}
И, это R код, который работает.
library(Rcpp)
library(RcppArmadillo)
library(RcppEigen)
library(RcppNumerical)
library(survival)
library(aftgee)
sourceCpp("C:/Users/mattw/Documents/paper_wj/exercise/aftsrr_wj_cpp/aftsrr_wj/optim_bfgs.cpp")
U_beta_r_non = function(beta,Time,Delta,Covari) {
n=length(Time)
p=ncol(Covari)
e_i_beta = as.vector(log(Time) + Covari %*% beta)
order_resid = order(e_i_beta)
Time = Time[order_resid]
Covari = matrix(Covari[order_resid,],nrow = n)
Delta = Delta[order_resid]
e_i_beta = e_i_beta[order_resid]
U_beta = list(NA)
for (i in 1:n) {
U_beta[[i]] = colSums(Delta[i]*t(Covari[i,]-t(Covari))*(e_i_beta>=e_i_beta[i]))
}
U_beta = Reduce('+',U_beta)/n
U_beta = sum(U_beta^2)
# grad = c(); h = 1e-4;
# for (it in 1:p) {
# beta_g = beta;
# beta_g[it] = beta[it] + h
#
# e_i_beta = as.vector(log(Time) + Covari %*% beta_g)
#
# order_resid = order(e_i_beta)
#
# Time = Time[order_resid]
# Covari = matrix(Covari[order_resid,],nrow = n)
# Delta = Delta[order_resid]
# e_i_beta = e_i_beta[order_resid]
#
# grad_beta = list(NA);
# for (i in 1:n) {
# grad_beta[[i]] = colSums(Delta[i]*t(Covari[i,]-t(Covari))*(e_i_beta>=e_i_beta[i]))
# }
# grad_beta = Reduce('+',grad_beta)/n
# grad_beta = sum(grad_beta^2)
#
# grad[it] = (grad_beta-U_beta)/h
# }
# return(U_beta)
return(sum(U_beta^2))
# return(grad)
}
#------------------------DATA GENERATION----------------------
set.seed(1)
n=300
beta_0=1
gamma_0=0.5
Z1=matrix(rnorm(n,3,1),nrow=n)
Z2=matrix(rexp(n,5),nrow=n)
T_aft=as.vector(exp(-beta_0*Z1-gamma_0*Z2+rnorm(n,5,1)))
C_aft=as.vector(exp(-beta_0*Z1-gamma_0*Z2+rnorm(n,6,1)))
X_aft=C_aft*(T_aft>C_aft)+T_aft*(T_aft<=C_aft)
D_aft=0*(T_aft>C_aft)+1*(T_aft<=C_aft)
table(D_aft)
beta_aftsrr=-aftsrr(Surv(X_aft,D_aft)~Z1+Z2)$beta;beta_aftsrr
init_beta = rep(0,2)
cpp_bfgs_esti = aftsrr_bfgs(init_beta,X_aft,D_aft,cbind(Z1,Z2),rep(0,2));cpp_bfgs_esti
optim_lbfgsb = optim(init_beta,function(x){U_beta_r_non(x,X_aft,D_aft,cbind(Z1,Z2))},method = "L-BFGS-B")$par;optim_lbfgsb
Теперь он работает, однако, кажется, что он просто дает init_beta без каких-либо вычислений в примере R.Поскольку некоторые люди говорят, что невозможно использовать вышеуказанные коды броненосцев, я изучаю библиотеку Эйгена.Тем не менее, я не знаком с вычислениями, трудно преобразовать их.Так что, надеюсь, есть способ использовать его с некоторыми изменениями.Спасибо!
Любые комментарии будут полезны.