I am fairly certain that this line is not kosher (and in any case it is a confusing):<div><br></div><div>  w = pd * 0.4 / (R::pnorm((qpd - sqrt(rsq) * (-0.42) * sgtemp), 0, 1, 1, 0));</div><div><br></div><div>try this instead:</div>
<div><br></div><div><div>w = qpd - sqrt(rsq) * (-0.42) * sgtemp</div><div>w = pd * 0.4 / (R::pnorm(w, 0, 1, 1, 0));</div><br><div class="gmail_quote">On Wed, Mar 13, 2013 at 11:42 PM, Aileen Lin <span dir="ltr"><<a href="mailto:aileenshanhong.lin@gmail.com" target="_blank">aileenshanhong.lin@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>My C code:</div>
<div>//[[Rcpp::depends("Rcpp")]]<br>#include <Rcpp.h><br>#include <iostream><br>using namespace Rcpp;</div>
<div> </div>
<div>//[[Rcpp::export]]<br>NumericVector sigmutest(double pd, double rsq){<br>  double qpd = R::qnorm(pd, 0, 1, 1, 0);<br>  double sgtemp = 0.2;<br>  double sg = 0.3;<br>  double eor = 1;<br>  <br>  double w = 0;<br>  while (eor>=0.0001) {  <br>

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

        std::cout << "sgtemp " << sgtemp << std::endl;<br>        eor = abs(sg - sgtemp);<br>        <br>        std::cout << "error " << eor << std::endl;<br>        <br>

      } <br>  NumericVector out(3);<br>  out(0) = sg;<br>  out(1) = sgtemp;<br>  out(2) = eor;<br>  return out;  <br>}<br clear="all"><br>My R code:</div>
<div><span style="border-collapse:separate;text-indent:0px;letter-spacing:normal;font:13px/15px 'Lucida Console';text-transform:none;white-space:pre-wrap;background-color:rgb(225,226,229);word-spacing:0px"><pre style="BORDER-BOTTOM-STYLE:none;LINE-HEIGHT:1.2;BORDER-RIGHT-STYLE:none;MARGIN:0px;OUTLINE-STYLE:none;FONT-FAMILY:'Lucida Console';WHITE-SPACE:pre-wrap!important;BORDER-TOP-STYLE:none;FONT-SIZE:10pt!important;BORDER-LEFT-STYLE:none">
<span style="COLOR:blue">> Rcpp::sourceCpp('src/sbi.cpp')
</span><span style="white-space:pre-wrap;COLOR:blue">> </span><span style="COLOR:blue">x <- sigmutest(0.0002327279, 0.1025499338)
</span>sg 0.2
sgtemp 0.219135
error 0</pre></span></div>
<div> </div>
<div>Does anyone know what is going on? Thanks.</div><span class="HOEnZb"><font color="#888888">
<div>-- <br>Aileen L.</div>
<div><br>
<div>View my Linkedin profile: <a style="padding-left:0px;padding-right:0px;border-bottom-width:0px;padding-top:0px;padding-bottom:0px;border-top-width:0px;vertical-align:baseline;line-height:15px;outline-style:none;color:rgb(0,102,153);font-size:13px;border-right-width:0px;margin:0px;font-family:Arial,Helvetica,'Nimbus Sans L',sans-serif;border-left-width:0px" title="View public profile" href="http://au.linkedin.com/in/aileen2" name="13d66fe227623061_SafeHtmlFilter_SafeHtmlFilter_webProfileURL" target="_blank">http://au.linkedin.com/in/aileen2</a></div>


<div><br></div>
<div>
<h2 style="padding-right:2px;padding-left:0px;padding-top:2px;text-align:left;font-family:Geneva,Arial,Helvetica,sans-serif;margin:0px;font-weight:normal;min-height:35px;padding-bottom:0px">
<font color="#993399">Being happy doesn't mean you're perfect. It just means you've decided to look beyond the imperfectio</font><font color="#cc33cc">ns- <a style="PADDING-BOTTOM:0px;MARGIN:0px;PADDING-LEFT:0px;PADDING-RIGHT:0px;TEXT-DECORATION:none;PADDING-TOP:0px" href="http://www.boardofwisdom.com/default.asp?topic=1010&search=K%2EB+Indiana+%28age+14%29" target="_blank">K.B Indiana (age 14)</a></font></h2>

</div></div>
</font></span><br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br></div>