multinomial on matrix row

[richfitz]

[richfitz]

ok, just worked out/remembered why this is hard - there's nothing that allows anything other than element wise updates on the left hand side in odin. So when we write

x[, ] <- y[i, j]

it's easy enough to imagine that this should become

for (i in 1..n) {
  for (j in 1..m) {
    x[i, j] <- y[i, y]
  }
}

but with a vector valued function this is hard! For a vector we allow

x[] <- rmultinom(n, p)

because that's unambiguous. but if we wanted to assign into a row or column of x there's no way of doing so:

x[, ] <- rmulitnom(n, p[, i])

could be updating x row-wise or column-wise. Getting this right probably requires more thought than I am likely to get to.

[richfitz]

Currently we allow

x[] <- rmultinom(n, p)

and we require that p and x have the same length. However, we should also allow things like:

m[, 1] <- rmultinom(n, p)

and require that dim(m, 1) == length(p) - this could be achieved by adding a "constraints" section to the IR, which I think could eventually remove some of the pain from the interpolation data processing and also pave the way for things like matrix multiplication (#38)

Doing this for a loop though is a bit hard, for the reasons indicated above. So if we wanted to update all of m with each row becoming a draw from the multinomial (i.e., as above) then we might do:

m[., ] <- rmultinom(n, p)

which would be equivalent to

m[, 1] <- rmultinom(n, p)
m[, 2] <- rmultinom(n, p)
...

for (i in 1..) {
  m[, i] <- rmultinom(n, p)
}

typically with a matrix loop we do:

m[, ] <- f(i, j)

but here the . would remove that index level from the loop so that to use different p and n values for each iteration we might write

m[., ] <- rmultinom(n[j], p[j])

(See anne's comment below - not quite right here as we need p to be a matrix...)

which would be equivalent to

m[, 1] <- rmultinom(n[1], p[1])
m[, 2] <- rmultinom(n[2], p[2])
...

for (j in 1..) {
  m[, j] <- rmultinom(n[j], p[j])
}

[richfitz]

For consistency, I think we should deprecate the current syntax for multinomials and require

x[.] <- rmultinom(n, p)

[annecori]

4 equation groups above, m[., ] <- rmultinom(n[j], p[j]) should be m[., ] <- rmultinom(n[j], p[,j])

[richfitz]

the other part of the puzzle is in #213