Ivan's answer is very impressive. Michael,
I am open to whatever suggestions but I am not aware a simple closed-form solution for my original question. Thanks, Zhu > On May 11, 2019, at 8:01 AM, Michael Weylandt <michael.weyla...@gmail.com> > wrote: > > This is very cool, but I wonder if it isn't over-kill for the larger > problem. > > In general, calculating the coefficient of an intercept-only GLM is just > calculating (a transformation of) the MLE of a univariate exponential > family distribution. (Things may be a bit trickier if the GLM also involves > weights and offsets, not just an intercept, but I'm assuming it doesn't.) > > OP: Can you clarify why you want to invoke R's entire GLM machinery as > opposed to just using the closed form solutions? > > Michael > >> On Sat, May 11, 2019 at 7:23 AM Ivan Krylov <krylov.r...@gmail.com> wrote: >> >> On Fri, 10 May 2019 16:17:42 +0000 >> "Wang, Zhu" <wan...@uthscsa.edu> wrote: >> >>> Are there any examples or links for me to follow through more closely? >> >> Calling R functions from C++ is described at >> <http://dirk.eddelbuettel.com/code/rcpp/Rcpp-quickref.pdf> and >> elsewhere in Rcpp documentation. An example follows: >> >> --------------8<--------------glmfit.cpp--------------8<-------------- >> #include <algorithm> >> #include <Rcpp.h> >> using namespace Rcpp; >> >> extern "C" double intercept_glm(size_t n, const double * response) { >> // access functions from default environment >> Function glm_fit("glm.fit"), coef("coef"); >> >> // intercept-only model: response ~ 1 >> NumericVector x(n); >> x.fill(1); >> >> // I couldn't find a way to wrap a double* into a NumericVector >> // without copying anything, sorry; perhaps someone else >> // can offer a solution >> NumericVector y(n); >> std::copy_n(response, n, y.begin()); >> >> // call the R function, convert the result back >> return as<double>(coef(glm_fit(x, y))); >> } >> --------------8<--------------glmfit.cpp--------------8<-------------- >> >> Since this function is extern "C" and uses only primitive C types, it >> should be fairly easy to call from Fortran. (C is the lingua franca of >> programming languages). Fortran-C interoperability is well described in >> "Modern Fortran Explained" by Metcalf et al. Here is the Fortran side >> of the code: >> >> --------------8<--------------callglm.f90--------------8<-------------- >> subroutine callglm(ret) >> use, intrinsic :: iso_c_binding, only: c_size_t, c_double >> ! using iso_c_binding here >> ! - to get correct type of ret when R calls the function >> ! - to convert variables before calling C function >> implicit none >> ! using F77-style arguments to match expectations of .Fortran() >> real(c_double), intent(out) :: ret >> ! toy data to compare against R code later >> real :: y(10) = [10, 11, 20, 9, 10, 8, 11, 45, 2, 3] >> ! the interface block declares an extern "C" function >> interface >> ! double intercept_glm(size_t n, const double * response) >> function intercept_glm(n, response) bind(c) >> use, intrinsic :: iso_c_binding >> real(c_double) :: intercept_glm >> integer(c_size_t), value :: n >> real(c_double) :: response(*) >> end function >> end interface >> >> ! call the function as you would call any other function >> ret = intercept_glm(int(size(y), c_size_t), real(y, c_double)) >> end subroutine >> --------------8<--------------callglm.f90--------------8<-------------- >> >> For a quick test, make sure that you have Rcpp installed and run: >> >> # adjust R version and path if your library is elsewhere >> PKG_CPPFLAGS='-g -I ~/R/x86_64-pc-linux-gnu-library/3.3/Rcpp/include' \ >> R CMD SHLIB callglm.f90 glmfit.cpp >> R >> library(Rcpp) >> dyn.load('callglm.so') # change extension if needed >> .Fortran('callglm', ret=numeric(1)) >> # $ret >> # [1] 12.9 >> coef(glm.fit(rep(1, 10), c(10, 11, 20, 9, 10, 8, 11, 45, 2, 3))) >> # [1] 12.9 >> >> To use this in a package, place both files in the src/ subdirectory of >> your package and add LinkingTo: Rcpp in the DESCRIPTION. >> >> -- >> Best regards, >> Ivan >> >> ______________________________________________ >> R-package-devel@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-package-devel >> > > [[alternative HTML version deleted]] > > ______________________________________________ > R-package-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-package-devel ______________________________________________ R-package-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-package-devel
