Reputation: 1775
What is the fastest way to do this summation in R?
This is what I have so far
ans = 0
for (i in 1:dimx[1]){
for (j in 1:dimx[2]){
ans = ans + ((x[i,j] - parameters$mu)^2)/(parameters$omega_2[i]*parameters$sigma_2[j])
}
}
where omega_2
, and sigma_2
are omega^2 and sigma^2 respectively.
Upvotes: 7
Views: 8756
Reputation: 48211
Usually it is quite easy to vectorize this kind of operations. In this case, though, it is a bit trickier when n
is not equal to m
and also because of double summation. But here is how we can proceed:
# n = 3, m = 2
xs <- cbind(1:3, 4:6)
omegas <- 1:3
sigmas <- 1:2
mu <- 3
sum((t((xs - mu) / omegas) / sigmas)^2)
# [1] 5
Here we use recycling three times and t()
to divide appropriate elements by sigmas
.
Upvotes: 3
Reputation: 44614
Nothing fancy:
# sample data
m <- matrix(1:20, 4)
sigma <- 1:ncol(m)
omega <- 1:nrow(m)
mu <- 2
sum(((m - mu) / outer(omega, sigma))^2)
Upvotes: 7