# R code is slow on filling a matrix

I am running an R code, which is implementing 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]))
p2=function(thetax)
dkde(thetax, fhatt)

fhatd=kde(x=deltat.mh[m1:m,s,t-1], h=hpi(deltat.mh[m1:m,s,t-1]))
p3=function(deltatx)
dkde(deltatx, fhatd)
# P1 and P2 are two distribution updated from the previous values of tetha.mh through kernel density estimation
h1=function(thetax)
log(p1(y[s,t],thetax)*p2(thetax))
# h1 is the the logarithm of the multiplication of p1 and p2 (laplace form estimate)
for(i1 in 2:5000)
{
thetas=rnorm(1,theta.mh[i1-1,s,t],stheta[i1,s,t])
r1[i1,s,t]=exp(h1(thetas)-h1(theta.mh[i1-1,s,t]))
if (r1[i1,s,t]>runif(1))
theta.mh[i1,s,t]=thetas
else {
theta.mh[i1,s,t]=theta.mh[i1-1,s,t]
}
}
-
What do you mean by "too much time"? –  ChrisW Jun 22 '12 at 0:13
It would help if you provide a reproducible example (code that I can paste directly into my console and run). –  bdemarest Jun 22 '12 at 0:35
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 –  Chase Jun 22 '12 at 1:12
How did you pre-allocate your matrix? –  Jason Morgan Jun 22 '12 at 3:01

## migrated from stackoverflow.comJun 22 '12 at 12:18

This question came from our site for professional and enthusiast programmers.

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...

-