Geek Answers Handbook

Matrix multiplication with Rcpp - different values ​​when assigning output

I wrote the following function to calculate the matrix (Weibull model information matrix)

#include <RcppEigen.h> #include <math.h> #include <vector> using namespace std; using Eigen::MatrixXd; // [[Rcpp::depends(RcppEigen)]] // [[Rcpp::export]] MatrixXd Weibull_FIM(const vector<double> x, const vector<double> w, const vector<double> param) { if(x.size() != w.size()){ Rcpp::Rcout<<"The length of x and w is not equal."<<std::endl; exit(1); } double a, b, lambda, h, constant; a = param[0]; b = param[1]; lambda = param[2]; h = param[3]; a = a + 0; //just to not get a warning Eigen::MatrixXd Fisher_mat(4, 4); size_t i; for(i=0; i < x.size(); i++) { constant = exp(-lambda * pow(x[i], h)); Eigen::MatrixXd f(4, 1); f(0, 0) = 1; f(1, 0) = -constant; f(2, 0) = b*pow(x[i], h)*constant; f(3, 0) = b*pow(x[i], h)*constant * lambda * log(x[i]); Fisher_mat = w[i] * f * f.transpose() + Fisher_mat; } return Fisher_mat; }

Now I want to calculate the matrix for some values ​​in R:

  Weibull_FIM(x=c(1, 1.239, 1.749371, 5), w = rep(.25, 4), param=c(1, 1, 1, 2)) # [,1] [,2] [,3] [,4] #[1,] 1.00000000 -0.157545687 0.210509365 0.037774174 #[2,] -0.15754569 0.045985587 -0.053326010 -0.004757122 #[3,] 0.21050936 -0.053326010 0.066320481 0.008736574 #[4,] 0.03777417 -0.004757122 0.008736574 0.002864707

This matrix is ​​fine. Now, if I first assign the output and then print it, I will have a different matrix (some of the elements will be changed to some large values!).

 res <- Weibull_FIM(x=c(1, 1.239, 1.749371, 5), w = rep(.25, 4), param=c(1, 1, 1, 2)) print(res) # [,1] [,2] [,3] [,4] #[1,] 1.00000000 4.727161e+180 0.210509365 2.267126e+161 #[2,] -0.15754569 5.104678e+199 -0.053326010 -4.757122e-03 #[3,] 0.21050936 2.079498e+64 0.066320481 8.736574e-03 #[4,] 0.03777417 -4.757122e-03 0.008736574 2.864707e-03

As you can see, the elements (1, 2), (1, 4), (2, 2) and (3, 2) get even larger values. I would be grateful for any help.

Update1 in response to @Ronald:

  • version.string R version 3.1.2 (2014-10-31)
  • Windows 7 x64
  • x86_64-w64-mingw32 platform
  • Rcpp_0.11.3, RcppEigen_0.3.2.3.0, tools_3.1.2
  • Rtools version 3.2.0.1948
+8
matrix r rcpp
source share
1 answer

You forgot to initialize Fisher_mat zeros before adding, so sometimes they can contain random garbage. Addition, for example, Fisher_mat.setZero(); before the for loop gives you consistent output.

Check out this link: Introduction to RcppEigen .

+9
source share

More articles:


All Articles