如何用R生成服从任意有限离散分布的随机数

我在网上找到了如何用matlab生成服从任意有限离散分布的随机数的方法，代码如下：
function out=drnd(p,n)
%out=drnd(p,n)
%copyright:rocwoods
%p is the probability distribution matrix;
%n is the number of the samples you want to generate;
%=================
a=cumsum(p(2,: ));
b=rand(n,1);
out=zeros(1,n);
for k=1:n
c=find(a<b(k));
if (isempty(c))
out(k)=p(1,1);
else
out(k)=p(1,c(end)+1);
end
end

举例说明:
比方说一个随机变量可以取1 2 3 4 5这5个值。每个值的概率分别为0.1 0.3 0.25 0.25 0.1
则P=[1 2 3 4 5;0.1 0.3 0.25 0.25 0.1]
生成1000个样本，如下:
out=drnd(P,1000)
%结果省略

当然在R中可以直接调用sample函数得到结果，但我尝试着编写Rcode编写相同的程序（因为个人觉得R和matlab有很多相似之处），
比如说X的分布律满足x=(x1,x2,.....,xk)，相应概率为p=(p1,p2,....,pk)，要生成n个满足所给分布律的随机数z，则先要生成n个（0,1）随机数，记为r，
当0<r[i]<=p1时，z[i]=x1; 当p1<r[i]<=p1+p2时，z[i]=x2；........  当pk<r[i]<=1时，z[i]=xk; 我编写了R程序，但是结果老是出错。
本人刚接触软件不久，希望高手赐教，悬赏不多，还请包涵。

drnd <- function(x, p, n) {
    z <- NULL
    ps <- cumsum(p)
    r <- runif(n)
    for (i in 1:n) z <- c(z, x[which(r[i] <= ps)[1]])
    return(z)
}
> x <- c(1, 2, 3, 4, 5)
> p <- c(0.1, 0.2, 0.2, 0.35, 0.15)
> out <- drnd(x, p, 10000)
> table(out)/10000
out
     1      2      3      4      5
0.0978 0.2056 0.1974 0.3480 0.1512

缺点：慢，因为那个for循环

等了两天了，你写的还不错，那分就给你吧，是不是我点最佳答案就可以了？
还有两个小问题，z <- c(z, x[which(r[i] <= ps)[1]])是什么意思，不是有两边吗，r[i]还要大于某某呀；
另外你最后搞个table(out)/10000是什么意思

哦，忘了说谢谢了

table(out)/10000是用来检查一下结果合理不
z <- c(z, x[which(r[i] <= ps)[1]])是找到离散的cdf里面第一个大于r[i]的数，然后z[i]就赋值为x里面这个位置的数。就是按照你帖子里的算法来的。

MUCHAS GRACIAS ！

De Nada  !

统计模拟，你可以参见武汉大学出版社《统计模拟及其R应用》，这本书挺好的

f<-function(n,x,prob){
  num<-length(prob)
  u<-runif(n)
  cup<-cumsum(prob)
  z<-x[as.integer(cut(u,breaks=c(0,cup)))]
  z}

一诺9257 发表于 2011-5-25 08:20
统计模拟，你可以参见武汉大学出版社《统计模拟及其R应用》，这本书挺好的

你说的是<<统计模拟与R实现>>吗?