矛盾した結果

Question

UPDATE 前の例では、したがって、以下に示すように、私は単純な例を使用できるようにしてください、複雑である：ここで

はRcppコードです：

#include <RcppArmadillo.h> 
#include <RcppArmadilloExtensions/sample.h> 
#include <Rmath.h> 
#include <Rcpp.h> 
// [[Rcpp::depends(RcppArmadillo)]] 
using namespace Rcpp ; 
using namespace arma; 
using namespace std; 

// [[Rcpp::export]] 
double chooseC(double n, double k) { 
    return Rf_choose(n, k); 
} 
// [[Rcpp::export]] 
double function3(double n, double m, double beta) { 
    double prob; 
    NumericVector k(m); 
    NumericVector k_vec(m); 
    if(n<m){prob=0;} 
    else{ 
    if(chooseC(n,m)==R_PosInf){ 
     k=seq_len(m)-1; 
     k_vec= (n-k)/(m-k)*std::pow((1-beta),(n-m)/m)*beta; 
      prob=std::accumulate(k_vec.begin(),k_vec.end(), 1, std::multiplies<double>())*beta; 
    } 
    else{ 
     prob = beta * chooseC(n,m) * std::pow(beta,m) * std::pow((1-beta),(n-m)); 
    } 

    } 
    return(prob); 
} 

ここ

はRコードです：

function4 <- function (n , m , beta) 
{ 
    if (n < m) 
    { 
    prob <- 0.0 
    } 
    else 
    { 
    if (is.infinite(choose(n,m))){ 
     k<-0:(m-1) 
     prob <- beta *prod((n-k)/(m-k)*(1-beta)^((n-m)/m)*beta) 
    } 
    else{ 
     prob <- beta * choose(n,m) * beta^m * (1-beta)^(n-m) 
    } 
    } 
    prob 
} 

比較：

input<-619 
beta<-0.09187495 

x<-seq(0, (input+1)/beta*3) 
yy<-sapply(x,function(n)function3(n,input, beta=beta)) 
yy2<-sapply(x,function(n)function4(n,input, beta=beta)) 
sum(yy)=0 
sum(yy2)=1 

しかし、他の入力を有する：

input<-1 
beta<-0.08214248 

両方の結果は、sum(yy)=sum(yy2)=0.9865887、同じです。

私はRcppとRコードの間の一貫性のない精度を引き起こす可能性が他に何かわからない、Rcppコードでdoubleを使用。

ありがとうございます！

Answer 1

私はRcppコードを修正したと思うので、結果が非​​常に小さい値の場合、RcppとRコードの両方で同じ結果が得られます。解決策は以下のように示されている：

#include <RcppArmadillo.h> 
#include <RcppArmadilloExtensions/sample.h> 
#include <Rmath.h> 
#include <Rcpp.h> 
// [[Rcpp::depends(RcppArmadillo)]] 

using namespace Rcpp ; 
using namespace arma; 
using namespace std; 

// [[Rcpp::export]] 
double chooseC(double n, double k) { 
    return Rf_choose(n, k); 
} 

// [[Rcpp::export]] 
double function3(double n, double m, double beta) { 
    double prob; 
    arma::vec k = arma::linspace<vec>(0, m-1, m); 
    arma::vec k_vec; 

    if(n<m){prob=0;} 
    else{ 
    if(chooseC(n,m)==R_PosInf){ 
     k_vec= (n-k)/(m-k)*pow((1-beta),(n-m)/m)*beta; 
     prob=arma::prod(k_vec)*beta; 
    } 
    else{ 
     prob = beta * chooseC(n,m) * pow(beta,m) * pow((1-beta),(n-m)); 
    } 

    } 
    return(prob); 
} 

このようにコードを記述することで高精度の矛盾を修正する理由しかし、私はまだ理解していません。 RcppとRcppArmadilloはまだ私には黒い箱のように見えます。