**主要问题:**如何在Rcpp中格式化一个复杂的结构体,使其能被optim.c接受?当返回到定义的optimfn函数时,如何修改此结构体?
**背景:**我正在Rcpp中编写一个EM算法,将用于beta混合模型。我试图通过R API(nmmin)使用一个简单的nelder mead数值优化方法,但是我很难通过变量列表。我对C++/Rcpp非常陌生。
信息:
eout计数:
- zprobs:矩阵(double[45203 x5])。(列和行的数量将变化)。
- parmlist:一个包含avec、mvec和svec的列表(每个列表有5个double)。(长度会有所不同)。
- xi:长度为45203的双精度数列表。(长度会有所不同)。
- denom:长度为45203的双精度数列表。(长度会有所不同)。
为了将我的变量传递给nmmin
,我尝试将它转换为结构体格式,首先将NumericMatrix和NumericVector转换为double *var
,(但我后来了解到这些不是真实的(编辑:这些不包含向量)),然后用这些项创建一个具有我定义的结构的对象。当前的方法能够编译,但我不确定*ex
进入enlikeB
函数时是什么样子。我还尝试将Rcpp::List
项转换为std::vectors
,以防这些项可能通过。
我现在在哪里:
#include <RcppCommon.h>
typedef struct{
double *xi;
int xil;
int zrows;
int zcols;
double *zdata;
double *denom;
int denoml;
}Ething;
#include <Rcpp.h>
//HELPER FUNCTIONS HERE:
//
// Defining the optimfn function:
typedef double optimfn(int n, double *par, void *ex);
double elnlikeB(int n, double *par, void *ex) {
Rcpp::Rcout << "entering elnlikeB\n";
Rcpp::Rcout << "The value of ex: " << ex << "\n";
Ething *et = (Ething *) ex;
Rcpp::Rcout << "The value of et : " << et << "\n";
Rcpp::Rcout << "making z\n";
int zrows = et ->zrows;
int zcols = et ->zcols;
double *z = et ->zdata;
Rcpp::NumericMatrix Zprobsmat;
int count = 0;
for(int x = 0; x < zcols; x++) {
for(int y = 0; y < zrows; y++) {
Zprobsmat(y, x) = z[count]; //row, column
count = count +1;
}
}
Rcpp::Rcout << "The value of zprob : " << Zprobsmat << "\n";
Rcpp::Rcout << "making xi\n";
int xil = et -> xil;
double *xid = et -> xi;
Rcpp::NumericVector xi;
for(int nn = 0; nn < xil; nn++){
xi[nn] = xid[nn];
}
Rcpp::Rcout << "The value of xi : " << xi << "\n";
Rcpp::Rcout << "making denom\n";
int denoml = et ->denoml;
double *denomd = et -> denom;
Rcpp::NumericVector denom;
for(int nn = 0; nn < denoml; nn++){
denom[nn] = denomd[nn];
}
Rcpp::Rcout << "The value of denom: " << denom << "\n";
Rcpp::Rcout << "making myguess\n";
Rcpp::NumericVector myguess;
for(int nn = 0; nn < n; nn++){
myguess[nn] = par[nn];
}
Rcpp::Rcout << "The value of myguess: " << myguess << "\n";
//This section works is the vectors and matrix are correct//
Rcpp::Rcout << "calculating\n";
int nnew = xi.size();
int nmixt = Zprobsmat.ncol();
int gsize = myguess.size();
Rcpp::NumericVector avec(nmixt);
avec[Rcpp::Range(0, (nmixt - 2))] = invmlogitc(myguess[Rcpp::Range(0, (nmixt - 2))]);
avec[(nmixt-1)] = (1 - rcpp_sum(avec));
Rcpp::NumericVector mvec(nmixt);
Rcpp::NumericVector svec(nmixt);
mvec = mutransform(myguess[Rcpp::Range((nmixt-1), (((nmixt-1)*2)))]);
svec = stransform(myguess[Rcpp::Range((((nmixt-1)*2)+1), (gsize-1))]);
Rcpp::NumericMatrix loglikemat( nnew, nmixt );
Rcpp::NumericMatrix lnfr( nnew, nmixt );
double sum = 0;
for(int i = 0; i < nmixt; i++){
double a = log(avec[i]);
double m = mvec[i];
double s = svec[i];
Rcpp::NumericVector ab = alphabetacalc(m, s);
for( int z = 0; z < nnew; z++){
lnfr(z,i) = (R::dbeta( xi[z], ab[0], ab[1], true));
double left = (Zprobsmat(z,i)/denom[z]);
double right = (a + lnfr(z,i));
loglikemat(z, i) = (left*right);
sum += loglikemat(z, i);
}
}
sum = (-sum);
Rprintf("sum: %d \n", sum);
return(sum);
}
optimfn elnlikeB;
//More helpers
extern "C" {
void nmmin(int n, double *xin, double *x, double *Fmin, optimfn fn,
int *fail, double abstol, double intol, void *ex,
double alpha, double beta, double gamma, int trace,
int *fncount, int maxit);
}
// [[Rcpp::export]]
Rcpp::List optim_test(Rcpp::List eout){
Rcpp::Rcout << "starting\n";
Rcpp::Rcout << "making avec, mvec, svec\n";
Rcpp::List parmlist = eout["parm.list"];
Rcpp::NumericVector avec = parmlist["avec"];
Rcpp::NumericVector mvec = parmlist["mvec"];
Rcpp::NumericVector svec = parmlist["svec"];
Rcpp::Rcout << "making zprob and xi \n";
Rcpp::NumericMatrix Zprobsmat = eout["zprob"];
Rcpp::NumericVector xi = eout["xi"];
int xil = xi.size();
Rcpp::NumericVector denom = eout["denom"];
int denoml = denom.size();
double nmixt = avec.size();
Rcpp::NumericVector mas = mlogitc(avec, nmixt);
Rcpp::Rcout << "making guess\n";
int gl = ((nmixt-1) + nmixt + nmixt);
Rcpp::NumericVector guess(gl);
guess[Rcpp::Range(0, (nmixt-2))] = mas;
guess[Rcpp::Range((nmixt-1), ((nmixt-1)*2))] = min(mvec);
guess[Rcpp::Range((((nmixt-1)*2)+1), (gl-1))] = sin(svec);
Rcpp::Rcout << "making zdata\n";
int tl = (xil*nmixt);
Rcpp::NumericVector zdata(tl);
int count = 0;
for(int x = 0; x < nmixt; x++) {
for(int y = 0; y < xil; y++) {
zdata[count] = Zprobsmat(y, x); //row, column
count = count +1;
}
}
int zrows = xil;
int zcols = nmixt;
Rcpp::Rcout << "making et\n";
Rcpp::List et = Rcpp::List::create(Rcpp::Named("xi") = xi, Rcpp::Named("xil") = xil,
Rcpp::Named("zrows") = zrows,
Rcpp::Named("zcols") = zcols,
Rcpp::Named("zdata") = zdata, Rcpp::Named("denom") = denom,
Rcpp::Named("denoml") =denoml);
Rcpp::Rcout << "making guess vec\n";
double vec[gl];
for(int nn = 0; nn < gl; nn++){
vec[nn] = guess[nn];
}
double opar[gl];
double Fmin = 0.0;
int fail = 0;
double abstol = 1.0e-8;
double intol = 1.0e-8;
double alpha = 1.0;
double beta = 0.5;
double gamma = 2.0;
int trace =0;
int fncount = 0;
int maxit = 500;
Rcpp::Rcout << "Attempting nmmin\n";
nmmin(gl, vec, opar, &Fmin,
elnlikeB, &fail, abstol, intol, &et, alpha, beta,
gamma, trace, &fncount, maxit);
Rcpp::List res = Rcpp::List::create(Rcpp::_["Fmin"] = Fmin, Rcpp::_["fail"]=fail);
return(res);`
}
1条答案
按热度按时间z5btuh9x1#
解决了!
我定义了一个类来代替
typedef struct
:在将数据转换为正确的格式后,我能够使
void *ex
对象等于:这允许您在创建对象时超越设置内存/大小,并以易于解释的方式传递格式化的对象。完整的代码将在几个月内在github上的
mgaynor1/nQuack
下提供。