ntr=nte=100
p=3

X.te = matrix(rnorm(nte*p),nrow=nte) # matrix of nte rows, p columns
f.te = apply (X.te^2,1,sum) # sum of squares of values
y.te = f.te + rnorm (nte, sd=3)

ks = c(1,5,10,20,50,100)
names (ks) = as.character(ks)

test.preds = array(data=NA,dim=c(nte,length(ks),200),dimnames=list(NULL, names(ks),NULL))
test.errs= test.preds

for (ind in 1:200){
    X.tr = matrix(rnorm(ntr*p),nrow=ntr) # matrix of ntr rows, p columns
    f.tr = apply (X.tr^2,1,sum) # sum of squares of values
    y.tr = f.tr + rnorm (ntr, sd=3)
    
    #pairwise distances   
    norm.te = apply (X.te^2, 1, sum) # apply function to rows of matrix
    norm.tr = apply (X.tr^2, 1, sum) 

    mat.norm.te = matrix (data=norm.te, nrow=nte, ncol=ntr, byrow=F) # each column is the norms of the test observations
    mat.norm.tr = matrix (data=norm.tr, nrow=nte, ncol=ntr, byrow=T) # each row is the norms of the test observations

    # matrix of distances
    dist.te.tr = mat.norm.te + mat.norm.tr - 2*X.te%*%t(X.tr) # matrix multiplication
    
    for (i in 1:nte) {
        neighbors = order (dist.te.tr[i,])
        for (k in 1:length(ks)){
            test.preds[i,k,ind] = mean(y.tr[neighbors[1:ks[k]]])
            test.errs[i,k,ind] = y.te[i]-test.preds[i,k,ind]
        }
    }
}    

# Examples of caluclations: 

#1. Mean of predictions for observation 3 with 5 neighbors
mean(test.preds[3,"5",])

#2. Mean of predictions for all test observations using 5 neighbors
apply(test.preds[,"5",],1,mean)

#3. Estimating the bias with 5 neighbors: 
apply(test.preds[,"5",],1,mean)-f.te