[R] Simulation
gb
gb at stat.umu.se
Wed Dec 6 23:08:26 CET 2000
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list