Q:

Storing a value in a matrix using a loop

My code is:
set.seed(100)
n <- 10
d <- rep(NA, n)
d[1] <- 0
y <- runif(n)
a <- 10

for (i in 2:length(y)) {
    d[i] <- d[i-1] + y[i-1]
} #This creates an interval with endpoints at my y uniform RVs

store.x <- NULL
for (j in 1:a) {
    x <- runif(1, min =0, max =sum(y))
    for (i in 1:length(y)) {
        if (x <= d[i+1] && x > d[i]) {
            store.x[j] <- i
            break
        }
    }
}  #This tells you which interval my x uniform RV is in

Now instead of store.x being a vector that tells me what intervals x falls into, I want it to be stored in the corresponding row and column of a matrix with the value of 1.  So for my first x since it falls into interval 7, my matrix will be all zeros except for a one in the first row, seventh column and a one in the seventh row, first column.
Any ideas on how to do this would be much appreciated!
Thanks

A:

You can more efficiently implement this algorithm using vectorized calls to runif(), cumsum(), and findInterval():
set.seed(1L); ## seed the PRNG for reproducible results
n <- 10L; ## length of x, y, and result matrix dimensions
y <- runif(n); ## produce random interval lengths
yb <- cumsum(c(0,y)); ## calculate interval boundaries
x <- runif(n,0,yb[length(yb)]); ## generate random x values
i <- findInterval(x,yb); ## find which intervals contain the x values
m <- matrix(0,n,n); ## init the result matrix
m[matrix(c(seq_len(n),i,i,seq_len(n)),ncol=2L)] <- 1; ## symmetric 1 assignment
m; ## print result
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    0    1    0    0    0    0    0    0     0
##  [2,]    0    0    1    0    0    0    0    0    0     0
##  [3,]    1    1    0    0    0    0    1    0    0     0
##  [4,]    0    0    0    1    0    0    0    0    1     0
##  [5,]    0    0    0    0    0    0    0    1    0     0
##  [6,]    0    0    0    0    0    1    0    0    0     0
##  [7,]    0    0    1    0    0    0    1    0    0     0
##  [8,]    0    0    0    0    1    0    0    0    0     1
##  [9,]    0    0    0    1    0    0    0    0    0     0
## [10,]    0    0    0    0    0    0    0    1    0     0

Technically, the way you've written the interval test x <= d[i+1] && x > d[i] means that you want the left boundary to be open and the right boundary to be closed for each interval. It appears that findInterval() has only recently added support for this variation on its logic, specifically in the r69814 development snapshot (which will eventually become R-3.3.0). So you probably don't have access to it yet, but when you do, you can pass left.open=T to get that behavior.

If we have a elements of x allowing for a != n then we have to decide how to map the index of each element in x to the index of the result matrix. The above solution assumes that a == n and that there is a direct mapping between the two index domains.
If we consider the elements of x to correspond to result matrix indexes 1:n in a cyclic fashion, then we can define xi to be the cyclic index, taking yi as i, i.e. the y interval index in which the x element landed. We can then aggregate all (xi,yi) hits to produce a count for each cell.
If you still want symmetry, then we can further accumulate the count for each (yi,xi), thus double-counting each hit, once for each of the two symmetric cells.
set.seed(1L); ## seed the PRNG for reproducible results
n <- 10L; ## length of y and result matrix dimensions
a <- 100L; ## length of x
y <- runif(n); ## produce random interval lengths
yb <- cumsum(c(0,y)); ## calculate interval boundaries
x <- runif(a,0,yb[length(yb)]); ## generate random x values
i <- findInterval(x,yb); ## find which intervals contain the x values
m <- matrix(0,n,n); ## init the result matrix
hit <- as.matrix(aggregate(n~xi+yi,cbind(xi=seq_len(n),yi=i,n=1),sum)); ## aggregate hits
m[hit[,c('xi','yi')]] <- m[hit[,c('xi','yi')]]+hit[,'n']; ## add n to xi,yi
m[hit[,c('yi','xi')]] <- m[hit[,c('yi','xi')]]+hit[,'n']; ## add n to yi,xi
m; ## print result
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    0    1    2    0    3    3    1    2     0
##  [2,]    0    2    2    2    1    1    3    3    1     0
##  [3,]    1    2    0    4    1    5    4    2    0     1
##  [4,]    2    2    4   10    1    2    1    1    3     2
##  [5,]    0    1    1    1    0    3    2    6    0     2
##  [6,]    3    1    5    2    3    2    3    5    1     0
##  [7,]    3    3    4    1    2    3    4    2    3     3
##  [8,]    1    3    2    1    6    5    2    2    3     2
##  [9,]    2    1    0    3    0    1    3    3    2     2
## [10,]    0    0    1    2    2    0    3    2    2     0

