Multinomial probabilities

If you need fit a probability using optim(), it is more secure to fit an inverse logit value and use a logit transformation after to ensure that the values stay always between 0 and 1:
L <- c(1.2)
p <- 1/(1+exp(L))
And then the probabilities of the two events are p and 1-p.

However it is more difficult in the case of multinomial distribution because the sum of the probabilities must be 1. Here is the correct way to do. For N events, you need N-1 parameters:
L <- c(10, -1, 3, 7, 0)
p <- 1/(1+exp(L))
p
[1] 4.539787e-05 7.310586e-01 4.742587e-02 9.110512e-04 5.000000e-01

The trick is to transform these probabilities into conditional probabilities:

p.multinomial <- function(p) c(p, 1)*c(1, sapply(seq_along(p), function(i) prod(1-p[1:i])))

p.multinomial(p)
[1] 4.539787e-05 7.310254e-01 1.275420e-02 2.333885e-04 1.279708e-01 1.279708e-01

sum(c(p, 1)*c(1, sapply(seq_along(p), function(i) prod(1-p[1:i]))))
[1] 1

To do the reverse:

inv.p.multinomial <- function(x) {
  p <- x[1]
  for (i in 2:(length(x)-1))
    p <- c(p, x[i]/(prod(1-p[1:(i-1)])))
  return(p)
}

x <- p.multinomial(p)
inv.p.multinomial(x)
[1] 4.539787e-05 7.310586e-01 4.742587e-02 9.110512e-04 5.000000e-01