R code is slow on filling a matrix

Question

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]
          }  
      }

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