[Rcpp-devel] Call R function

Douglas Bates bates at stat.wisc.edu
Thu Jun 21 14:54:41 CEST 2012


On Wed, Jun 20, 2012 at 4:30 PM, bbonit at tin.it <bbonit at tin.it> wrote:
>
> i must use it to evaluate  matropolis ratio (green ratio) in the following
> code in R :  i have  start just  now to study c and c++   (like you can see
> by my stupids questions (for you surely))  if  you have some suggestion im
> very happy to "heard" all suggestion ......thank for all  of the mailing
> list .......that teach me much
> tahnk you again and again
> OneTimeVariable<-function (logpost, start, n.iter, burn, thin, scale, ...)
> {
>     pb <- txtProgressBar(min = 0, max = n.iter, style = 3)
>     p = length(start)
>     vth = array(0, dim = c(n.iter/thin, p))
>     f0 = logpost(start, ...)
>     arate = array(0, dim = c(1, p))
>     th0 = start
>     for (i in (-burn):n.iter) {
>         setTxtProgressBar(pb, i)
>         for (j in 1:p) {
>             th1 = th0
>             th1[j] = th0[j] + rnorm(1) * scale[j]
>             f1 = logpost(th1, ...)
>             u = runif(1) <= min(exp(f1 - f0),1)
>             th0[j] = th1[j] * (u == 1) + th0[j] * (u == 0)
>             f0 = f1 * (u == 1) + f0 * (u == 0)
>             vth[floor(i/thin), j] = th0[j]
>             arate[j] = arate[j] + u
>         }
>     }
>     arate = arate/n.iter
>     stuff = list(par = mcmc(vth), accept = arate)
>     return(stuff)
> }

In the enclosed file I provide two alternative versions in R and one
using Rcpp but none of these are tested.  I suggest that you create a
test case and first find out if they produce the same distribution of
results (OneTimeV1 should be faster but the results will be different,
even after calling set.seed with the same seed).  Then use the
benchmark package to try see the relative speed.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: foo.R
Type: application/octet-stream
Size: 4243 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20120621/42bf0363/attachment.obj>


More information about the Rcpp-devel mailing list