矩阵乘法 - 不同的值时，输出分配

我写了下面的函数来计算的矩阵（维泊尔模型的信息矩阵）矩阵乘法 - 不同的值时，输出分配

#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; 
}

现在我想要计算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

该矩阵是正确的。现在，如果我首先分配输出然后打印它，我会得到一个不同的矩阵（一些元素将被改为一些大的值！）。

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

如可以看到的，元件（1,2），（1,4），（2,2）和（3,2）得到另一个大的值。我将不胜感激任何帮助。

UPDATE1响应于@Ronald：

version.stringř版本3.1.2（2014年10月31日）
视窗7 64
平台x86_64的-W64-的mingw32
Rcpp_0.11.3，RcppEigen_0.3.2.3.0，tools_3.1.2
Rtools版本3.2.0.1948

来源

2015-04-22 Ehsan Masoudi

您应该提供您的系统信息。我可以在Win7系统上使用R 3.2.0以及最新版本的Rtools（我需要我们的系统管理员进行更新）来重现这一点。我无法在具有R 3.1.2和最新软件包版本的Linux系统上重现此操作。 – Roland

在总结之前忘记初始化Fisher_mat零，因此它有时可能包含随机垃圾。添加例如Fisher_mat.setZero();之前for循环会给你一个一致的输出。

查看此参考文献：RcppEigen Introduction。

来源

2015-04-22 12:56:06 tonytonov

添加我建议的结果，两次结果都是正确的（我在一个非常类似的设置上检查过）。 – tonytonov

谢谢你，现在的结果是一致的。但如果我不设置它，那么Fisher_mat包含非常小的数字（如1.780215e-306），在R 1.780215e-306 + .001 = .001和1.780215e-306 - 中。 001 = -.001。好吧，让我们假设有精确的问题。但是我的主要问题是为什么（Weibull_FIM（x = c（1,1.239,1.749371,5），w = rep（.25,4），param = c（1,1,1,2）））给了我右边的答案但这一个不是（res < - Weibull_FIM（x = c（1,1.239,1.749371,5），w = rep（.25,4），param = c（1,1,1,2）））你的回答解决了我的问题（一个），但我一直在等待原因（在添加你的线路之前）。 –

@EhsanMasoudi它可能包含这些数字在特定情况下，但这些情况不保证是可重复的。这当然与精度无关。 – Roland

矩阵乘法 - 不同的值时，输出分配

回答

相关问题