Compute a double sum in R

Question

I have to compute a test statistic with a double sum.

I solved it like this:

T<-numeric(1)
for(j in 1:n){
   for(k in 1:n){
 T = T + ((1/n)*(exp(-(1/2)*((Y[j]-Y[k])^2))))}
 T = T - ((sqrt(2))*(exp(-(1/4)*((Y[j])^2))))}
 T = T + (n*(3^(-(1/2))))

Is there an easier way to compute the test statistic?

Answer 1

It's more useful to create the indexes in advance and then just sum over an array rather than computing new indices over two nested loops

indexes = expand.grid(1:n,1:n)
T = 1/n*sum(exp(-1/2*(Y[indexes[,1]]-Y[indexes[,2]])))
T = T-(sqrt(2))*sum(exp(-1/4*(Y[1:n])))
T = T+n/sqrt(3)

Edit: For large n, this is impractical, as an n of 1,000,000 would make a 3.7 TB data frame with expand.grid. You can always use the for loops, even if they are slow, but I would recommend using C++ if you need to have absurdly large N, because that is 1 trillion loops, which will take a very long time to compute.

Answer 2

Use

n=100;
Y=runif(100);
T=0;

Ydiff=outer(Y,Y,"-")^2;
Y_1=exp(-0.5*Ydiff);
Y_2=sqrt(2)*exp(-0.25*Y^2);

T=sum(rowMeans(Y_1)-Y_2) + (n*(3^(-(1/2))))

Comparison of methods given so far give:

T=0;
n=100;
set.seed(100)
Y=runif(100);
for(j in 1:n){
  for(k in 1:n){
    T = T + ((1/n)*(exp(-(1/2)*((Y[j]-Y[k])^2))));
  }
  T = T - ((sqrt(2))*(exp(-(1/4)*((Y[j])^2))));
}
T = T + (n*(3^(-(1/2))));
print(T)
#21.18983

T=0;
Ydiff=outer(Y,Y,"-")^2;
Y_1=exp(-0.5*Ydiff);
Y_2=sqrt(2)*exp(-0.25*Y^2);
T=sum(rowMeans(Y_1)-Y_2) + (n*(3^(-(1/2))));
print(T)
# 21.18983

T=0;
indexes = expand.grid(1:n,1:n);
T = 1/n*sum(exp(-1/2)*((Y[indexes[,1]]-Y[indexes[,2]])));
T = T-(sqrt(2))*sum(exp(-1/4*(Y[1:n])));
T = T+n/sqrt(3);
print(T)
# -66.71403