[R-pkg-devel] How to obtain intercept of intercept-only glm in Fortran?
Wang, Zhu
w@ngz1 @end|ng |rom uth@c@@@edu
Sat May 11 15:28:41 CEST 2019
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.weylandt using 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.r00t using gmail.com> wrote:
>>
>> On Fri, 10 May 2019 16:17:42 +0000
>> "Wang, Zhu" <wangz1 using 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 using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-package-devel
>>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-package-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-package-devel
More information about the R-package-devel
mailing list