From: Michael Braun <braunm_at_MIT.EDU>

Date: Thu 04 Jan 2007 - 21:09:36 GMT

Michael Braun

Assistant Professor of Marketing

MIT Sloan School of Management

38 Memorial Drive, E56-329

Cambridge, MA 02139

braunm@mit.edu

(617) 253-3436

R-devel@r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri Jan 05 08:12:17 2007

Date: Thu 04 Jan 2007 - 21:09:36 GMT

I am experiencing some odd behavior with the .Call interface, and I am hoping someone out there can help me with it. Below is a simple example (in that there are R packages that do exactly what I want), but this code illustrates the problem. Thanks in advance for any help you can provide.

x<-c(0,0,0,0,0,0) mu<-c(0,0,0,0,0,0) sig<-diag(6) print(sig) w<-.Call("R_mvnorm_logpdf",as.double(x),as.double(mu),sig, as.integer(6)) print(sig) # sig has changed after .Call

This code takes the SEXP's that were passed from R and converts them to gsl objects, and then calls the logpdf function:

# include <R.h> # include <Rinternals.h> # include <Rmath.h> # include <gsl/gsl_matrix.h> # include <gsl/gsl_vector.h> # include <gsl/gsl_blas.h> # include <gsl/gsl_linalg.h> # include <gsl/gsl_math.h> SEXP R_mvnorm_logpdf (SEXP xx, SEXP mux, SEXP sigmax, SEXP kx) { int k = INTEGER(kx)[0]; double * xAr = REAL(xx); double * muAr = REAL(mux); double * sigmaAr = REAL(sigmax); SEXP res; gsl_vector_view xView = gsl_vector_view_array(xAr,k); gsl_vector_view muView = gsl_vector_view_array(muAr,k); gsl_matrix_view sigmaView = gsl_matrix_view_array(sigmaAr,k,k); gsl_vector * x = &xView.vector; gsl_vector * mu = &muView.vector; gsl_matrix * sigma = &sigmaView.matrix; 1: double logans = gsl_MB_mvnorm_logpdf(x, mu, sigma, k); // <-call logpdf here PROTECT(res=allocVector(REALSXP,1)); REAL(res)[0] = logans; UNPROTECT(1); return(res); }

The logpdf function is here

double gsl_MB_mvnorm_logpdf(gsl_vector * beta, gsl_vector * betaMean, gsl_matrix * sigma, int k) { // computes density of multivariate normal vector at vector beta, with mean betaMean and cov sigma double logdetSigma = 0; double res; double * kern; int i, err; // pointer to Cholesky decomp of sigma gsl_matrix * sigmaChol = gsl_matrix_alloc(k, k); // define matrix that will store Chol decomp gsl_matrix_memcpy(sigmaChol, sigma); gsl_linalg_cholesky_decomp(sigmaChol); // compute logdet of sigma by 2*sum of log of diagomal elements of chol decomp for (i=0; i<k; i++) { logdetSigma = logdetSigma + log(gsl_matrix_get(sigmaChol,i,i)); } logdetSigma = 2*logdetSigma; // compute (beta-mean)' sigma^(-1) (beta-mean) gsl_vector * x = gsl_vector_alloc(k);

2: // gsl_matrix_fprintf(stdout,sigma,"%f");

gsl_vector_memcpy(x, beta); gsl_vector_sub(x, betaMean); // beta - betaMean gsl_vector * y = gsl_vector_alloc(k); gsl_vector_memcpy(y,x); gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, sigmaChol, y); // y = inv(chol)*x from BLAS gsl_blas_ddot(y,y,kern); // kern = y'y // compute log density res = -k*M_LN_SQRT_2PI - 0.5*(logdetSigma + *kern); // release space gsl_matrix_free(sigmaChol); gsl_vector_free(x); gsl_vector_free(y); return(res); } // end gsl_mvnorm_pdf

The problem is that after I make the .Call in R, the value of the sig matrix changes (the 1 in the upper left corner changes from a 1 to a 0). Since I don't make changes to the sigma object directly, I don't know how that could happen. The following output is the sig matrix, before and after the .Call.

> source("pdfTest.R") [,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 0 0 0 0 0[2,] 0 1 0 0 0 0[3,] 0 0 1 0 0 0[4,] 0 0 0 1 0 0

[5,] 0 0 0 0 1 0

[6,] 0 0 0 0 0 1

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 0 0 0 0 0 0[2,] 0 1 0 0 0 0[3,] 0 0 1 0 0 0[4,] 0 0 0 1 0 0

[5,] 0 0 0 0 1 0

[6,] 0 0 0 0 0 1

Through some debugging, I do know that it occurs somewhere in the logpdf function (line 1: in the first func). I should also note that this function does compute the density correctly.

Here's some other, potentially relevant information. Note the print statement in line 2 of the logpdf function. If that print were "uncommented" and placed somewhere *earlier* in the function, the elements of the sigma matrix are as they should be, and the program runs with no problem (except for the change in sig as described above). However, if the print statement is at line 2: or later, the correct matrix elements are printed, but then the .Call crashes with the following message:

- caught segfault *** address 0x6, cause 'memory not mapped'

(If I change the size of the matrix to diag(k), the 0x6 becomes 0xk). So, for some reason, k is interpreted as a memory address. I double-checked for places where I may have confused pointers and values, but I can't see where that would have happened. Also, after searching through the R-devel archives, I noticed that others have had other odd 'memory not mapped' errors in totally different scenarios.

So I am flummoxed, and don't really know where to go from here.

Best wishes,

MB

Michael Braun

Assistant Professor of Marketing

MIT Sloan School of Management

38 Memorial Drive, E56-329

Cambridge, MA 02139

braunm@mit.edu

(617) 253-3436

R-devel@r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri Jan 05 08:12:17 2007

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Thu 04 Jan 2007 - 22:31:01 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*