user1871528
user1871528

Reputation: 1775

Fastest way to do this double summation?

What is the fastest way to do this summationenter image description here 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

Answers (2)

Julius Vainora
Julius Vainora

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

Matthew Plourde
Matthew Plourde

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

Related Questions