Question: I would like to compute in R something of the following kind:

![formular](https://chart.googleapis.com/chart?cht=tx&chl=%5Csum_%7Bx%3D1%7D%5En%20%5Csum_%7By%3D0%7D%5Ex%20f(y,x)) 

It is important that the summand could be any function f(y,x).

My approach so far is with nested for loops:

```
n <- 5
fun <- function(y,x){y^2 + sqrt(y*x)} # might be any function of y and x
sum_x <- c()
for(x in 1:n){
  sum_y <- c()
  for(y in 0:x){
    sum_y[y+1] <- fun(y,x)
  }
  sum_x[x] <- sum(sum_y)
}
sum(sum_x) # 147.6317

```

However, I do not like this approach. It's pretty ugly and becomes very inconvenient if lower and upper bound need to be more flexible. I thought about using `expand.grid` and then applying `fun` to it using `mapply`, but couldn't figure out how to express the nested structure of the sums. Any suggestions how to do this?

Answer: 
You can perform the `outer` product based on a function. This outer product will look at all combinations of two input variables and place the result in a matrix; it takes the following form:

```
outer(<rows>, <cols>, FUN)

```

In your case specifically, the following suffices:

```
n <- 5

fun <- function(x, y) {ifelse (y > x, 0, y^2 + sqrt(x * y))}

outer(1:n, 1:n, FUN = fun) %>% sum() # 147.6317

```

Since `y` ranges from `0` and `y` occurs in both terms, it defaults to `0` (by chance). Regardless, it's necessary in this case to include some form of indexing in the function definition since the nested summing of `y` is dependent on `x`.
We can also use `outer` with use of `rowCumsums` from `matrixStats`

```
library(matrixStats)
sum(outer(seq_len(n), seq_len(n), FUN = fun) * rowCumsums(diag(n)))
#[1] 147.6317

```

---

Or with `crossing` from `tidyr`

```
library(tidyr)
library(dplyr)
crossing(x = seq_len(n), y = seq_len(n)) %>%
   filter(y <= x) %>%
   transmute(out = fun(y, x)) %>%
   summarise(out = sum(out)) %>% 
   pull(out)
#[1] 147.6317

```