[Rcpp-devel] not sure if it's a rcpp question or a cpp question
Alon Honig
honeyoak at gmail.com
Thu Mar 14 05:00:02 CET 2013
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:
> 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 at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130314/4150d86b/attachment.html>
More information about the Rcpp-devel
mailing list