[R] Simulation

[gb]

I got answers from Peter D., Ben Bolker, and Yudi Pawitan:

Peter's suggestion:
tabulate(sample(n, k, replace = T), n)
'tabulate' was what I couldn't find!

Yudi's is in the same spirit:
table(ceiling(n*runif(k)))
Change that to 'tabulate(ceiling(n*runif(k)), n)'
and I also get what I want.

Ben's
rmulti<-function (n, size, prob)
[............] (see his letter)

is much slower for _one_ replicate (n=1), but it vectorizes
and performs very well for large  n (note that Peter's and
Yudi's  'n'  is  'length(prob)'  in Ben's function).
And it puts no restriction on 'prob'.

We can however generalize  'sample' to
tabulate(sample(n, k, replace = T, prob = p), n)
but the vectorization is still to be done.

Many thanks to the three of you!

Göran

Original question:

> I want to draw a random sample Y_1, ..., Y_n, iid Poisson,
> but condition on their sum being equal to  k. Two
> suggestions:
>
> 1)   lambda <- k / n ## optimal value
>      ysum <- -1
>      while (ysum != k){
>         y <- rpois(n, lambda)
>         ysum <- sum(y)
>      }
>      y
>
> 2) Aiming at Multinom(k, 1/n, ..., 1/n):
>
>    y <- as.vector(table(factor(sample(n, k, replace = T), levels = 1:n)))
>
> 1) seems to be faster than 2), but my question is whether there is
> something even better?  'sample' is fast enough, but how can I do
> the rest (factor, table) smarter in 2)?

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._