2015-04-22 3 views
8

ho scritto la seguente funzione per calcolare una matrice (matrice di informazione di modello di Weibull)moltiplicazione di matrici con Rcpp - valori diversi quando l'uscita è assegnata

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

Ora voglio calcolare la matrice per alcuni valori 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 

Questa matrice va bene. Ora, se prima assegno l'output e poi lo stampo, avrò una matrice diversa (alcuni degli elementi saranno cambiati in alcuni grandi valori!).

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 

Come si vede, gli elementi (1, 2), (1, 4), (2, 2) e (3, 2) ottenere un altro grandi valori. Apprezzerei qualsiasi aiuto.

Update1 in risposta a @Ronald:

  • version.string R versione 3.1.2 (2014-10-31)
  • Windows 7 x64
  • piattaforma x86_64-W64-mingw32
  • Rcpp_0.11.3, RcppEigen_0.3.2.3.0, tools_3.1.2 versione
  • Rtools 3.2.0.1948
+1

È necessario fornire le informazioni di sistema. Posso riprodurlo usando R 3.2.0 con le versioni più recenti del pacchetto su un sistema Win7 con una versione precedente di Rtools (ho bisogno del nostro amministratore di sistema per l'aggiornamento). Non riesco a riprodurlo su un sistema Linux con R 3.1.2 e le più recenti versioni del pacchetto. – Roland

risposta

9

Hai dimenticato di inizializzare lo Fisher_mat con gli zeri prima della sommatoria, quindi a volte può contenere una spazzatura casuale. Aggiungendo per es. Fisher_mat.setZero(); prima di un ciclo for fornirà un risultato coerente.

Controllare questo riferimento: RcppEigen Introduction.

+0

Aggiungere la riga che propongo e il risultato sarà corretto entrambe le volte (l'ho verificato su una configurazione molto simile). – tonytonov

+0

Grazie, il risultato è coerente ora. Ma se non lo metto, il Fisher_mat contiene numeri molto piccoli (come 1.780215e-306) e in R 1.780215e-306 + .001 = .001 e 1.780215e-306 -. 001 = -.001. Ok supponiamo che ci sia stato un problema di precisione. Ma la mia domanda principale è perché (Weibull_FIM (x = c (1, 1.239, 1.749371, 5), w = rep (.25, 4), param = c (1, 1, 1, 2))) mi dà il diritto rispondi ma questo non è (res <- Weibull_FIM (x = c (1, 1.239, 1.749371, 5), w = rep (.25, 4), param = c (1, 1, 1, 2))) la tua risposta ha risolto il mio problema (uno su) ma continuo ad aspettare il motivo (prima di aggiungere la tua linea). –

+2

@EhsanMasoudi Potrebbe contenere questi numeri in circostanze specifiche, ma queste circostanze non sono garantite per essere riproducibili. Non è certamente correlato alla precisione. – Roland