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....@gmail.com>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. > -- > 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)<http://www.boardofwisdom.com/default.asp?topic=1010&search=K%2EB+Indiana+%28age+14%29> > > _______________________________________________ > Rcpp-devel mailing list > Rcpp-devel@lists.r-forge.r-project.org > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel >
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel