I am running a Markov Chain Monte Carlo algorithm for updating a density distribution. There is a specific section of my code which tries to fill a very large matrix thetha.mh (of dimension 5000x2x60). the following code snippet is trying to fill the first dimension of my matrix, but the specific operation in second line of the loop is taking too much time.

Does anybody know of a better way to improve the efficiency of this operation?

fhatt=kde(x=theta.mh[m1:m,s,t-1], h=hpi(theta.mh[m1:m,s,t-1]))
  dkde(thetax, fhatt)

fhatd=kde(x=deltat.mh[m1:m,s,t-1], h=hpi(deltat.mh[m1:m,s,t-1]))
  dkde(deltatx, fhatd)
# P1 and P2 are two distribution updated from the previous values of tetha.mh through kernel density estimation
# h1 is the the logarithm of the multiplication of p1 and p2 (laplace form estimate)
for(i1 in 2:5000)
        if (r1[i1,s,t]>runif(1)) 
        else {
  • \$\begingroup\$ What do you mean by "too much time"? \$\endgroup\$
    – ChrisW
    Jun 22, 2012 at 0:13
  • \$\begingroup\$ It would help if you provide a reproducible example (code that I can paste directly into my console and run). \$\endgroup\$
    – bdemarest
    Jun 22, 2012 at 0:35
  • \$\begingroup\$ other generic comments include compiling your function, moving anything out of the loop that you can, and consider using Rcpp and inline to rewrite in C++ depending on your wizard skills \$\endgroup\$
    – Chase
    Jun 22, 2012 at 1:12
  • \$\begingroup\$ How did you pre-allocate your matrix? \$\endgroup\$
    – Jason Morgan
    Jun 22, 2012 at 3:01

1 Answer 1


The functions you use (kde, dkde, hpi, ...) do not seem to be standard R functions and you are not very explicit on your exact data structur, so it's hard to help. Could you provide more info on that (packages used, ...)?

What I think to realize though - but correct me if I'm wrong - is the following. When substituting h1 in the second line of your loop with its function body, you could simplify the function and get rid of all exp and log, which should speed up your script at least a bit:

    exp( h1(thetas) - h1(theta.mh[i1-1,s,t]) )
--> exp( log(p1(...)*p2(...)) - log(p1(...)*p2(...)))
--> exp(log(p1(...)*p2(...))) / exp(log(p1(...)*p2(...)))
--> p1(...)*p2(...) / p1(...)*p2(...)

that means, what your second line actually computes is:

p1(y[s,t],thetas)*p2(thetas) / p1(y[s,t],theta.mh[i1-1,s,t])*p2(theta.mh[i1-1,s,t])

I hope that is of some use to you...


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy