# set-up
n = 100; m = 50; p = 5
gr = matrix(stats::rnorm(m*p), m, p)
L = R = matrix(stats::rnorm(n*p), n, p)
s1 = matrix(1, n, p)
missing.idx = sample.int(n = n*p, size = p*p)
L[missing.idx] = L[missing.idx] - stats::runif(p, 0, 1)
# R solution
lik = matrix(nrow = n, ncol = m)
for (i in 1:n) {
for(k in 1:m) {
lik[i,k] = prod(ifelse(
L[i,] == R[i,],
stats::dnorm(L[i,]-gr[k,], sd = s1[i,]),
stats::pnorm(R[i,]-gr[k,], sd = s1[i,]) - stats::pnorm(L[i,]-gr[k,], sd = s1[i,])
))
}
}
# Compare R to RcppParallel method
all.equal(lik, likMat(L, R, gr, s1))
Run the code above in your browser using DataLab