Skip to content

[Rcpp-devel] not sure if it's a rcpp question or a cpp question

3 messages · Aileen Lin, Dirk Eddelbuettel, Alon Honig

Original
Aileen Lin
# 
My C code:
//[[Rcpp::depends("Rcpp")]]
#include <Rcpp.h>
#include <iostream>
using namespace Rcpp;

//[[Rcpp::export]]
NumericVector sigmutest(double pd, double rsq){
  double qpd = R::qnorm(pd, 0, 1, 1, 0);
  double sgtemp = 0.2;
  double sg = 0.3;
  double eor = 1;

  double w = 0;
  while (eor>=0.0001) {
        sg = sgtemp;
        w = pd * 0.4 / (R::pnorm((qpd - sqrt(rsq) * (-0.42) * sgtemp), 0,
1, 1, 0));
        sgtemp = (-0.5) * w + 0.4;
        std::cout << "sg " << sg << std::endl;
        std::cout << "sgtemp " << sgtemp << std::endl;
        eor = abs(sg - sgtemp);

        std::cout << "error " << eor << std::endl;

      }
  NumericVector out(3);
  out(0) = sg;
  out(1) = sgtemp;
  out(2) = eor;
  return out;
}

My R code:
sgtemp 0.219135
error 0


Does anyone know what is going on? Thanks.
Dirk Eddelbuettel
#
On 14 March 2013 at 14:42, Aileen Lin wrote:
| My C code:
| //[[Rcpp::depends("Rcpp")]]
| #include <Rcpp.h>
| #include <iostream>
| using namespace Rcpp;
| ?
| //[[Rcpp::export]]
| NumericVector sigmutest(double pd, double rsq){
| ? double qpd = R::qnorm(pd, 0, 1, 1, 0);
| ? double sgtemp = 0.2;
| ? double sg = 0.3;
| ? double eor = 1;
| ?
| ? double w = 0;
| ? while (eor>=0.0001) {?
| ??????? sg = sgtemp;
| ??????? w = pd * 0.4 / (R::pnorm((qpd - sqrt(rsq) * (-0.42) * sgtemp),
| 0, 1, 1, 0));
| ??????? sgtemp = (-0.5) * w + 0.4;
| ??????? std::cout << "sg " << sg << std::endl;
| ??????? std::cout << "sgtemp " << sgtemp << std::endl;
| ??????? eor = abs(sg - sgtemp);
| ???????
| ??????? std::cout << "error " << eor << std::endl;
| ???????
| ????? }
| ? NumericVector out(3);
| ? out(0) = sg;
| ? out(1) = sgtemp;
| ? out(2) = eor;
| ? return out;?
| }
| 
| My R code:
| 
| > Rcpp::sourceCpp('src/sbi.cpp')
| > x <- sigmutest(0.0002327279, 0.1025499338)
| sg 0.2
| sgtemp 0.219135
| error 0
| 
| ?
| Does anyone know what is going on? Thanks.

abs() is a C function for int-on-int. So eor becomes 0. 

Use std::fabs() instead.

Dirk

| --
| Aileen L.
| 
| View my Linkedin profile:?http://au.linkedin.com/in/aileen2
| 
| 
| Being happy doesn't mean you're perfect. It just means you've decided to look
| beyond the imperfections-?K.B Indiana (age 14)
| 
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Alon Honig
# 
I am fairly certain that this line is not kosher (and in any case it is a
confusing):

  w = pd * 0.4 / (R::pnorm((qpd - sqrt(rsq) * (-0.42) * sgtemp), 0, 1, 1,
0));

try this instead:

w = qpd - sqrt(rsq) * (-0.42) * sgtemp
w = pd * 0.4 / (R::pnorm(w, 0, 1, 1, 0));

On Wed, Mar 13, 2013 at 11:42 PM, Aileen Lin
<aileenshanhong.lin at gmail.com>wrote:

            

      
      
      
    

        

      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130314/4150d86b/attachment.html>