本文在Creative Commons许可证下发布
思路如下：
首先需要设定一个红包金额分配机制。这里使用一个简化模型，即假设抢到红包的人分得的金额占总金额的比例（每个人对应一个比例，在这里设定为向量，且其总和为1）服从Dirichlet分布，参数为（α，α，…, α），α决定红包发放的“公平”程度，α越大，每人分得的金额比例就越倾向于平均，反之，波动性越大。好，下面是模拟程序：
设定红包接力100次，计算大家最后手上的余额。
library("gtools", lib.loc="D:/Program Files/R/R-3.3.1/library")
param=list(group_size=50,##群的人数
init_balance=50,##每人初始金额
hb_amount=20,##每次红包总金额
hb_size=10,##每次红包发多少人
niter=100,##红包接力次数
alpha=2##红包金额分配参数
)
hb_experiment=function(p)
{
id=1:p$group_size #分配群成员ID编号
balance=rep(p$init_balance,p$group_size)#当前每人资产
bal=matrix(0,p$niter,p$group_size)
for(i in 1:p$niter)
{
#破产就别玩了
players=id[balance>0]
if(i == 1)
{
#第一轮随机挑人发红包
host=sample(players,1)
}else{
#后面就找抢的最多的
host=winners[which.max(winner_amount)]
}
#红包主掏钱
balance[host]=balance[host]-p$hb_amount
#手快有，手慢无
winners=sample(players,p$hb_size)
#每人领取的红包金额
winner_amount=p$hb_amount*c(rdirichlet(1,rep(p$alpha,p$hb_size)))
balance[winners]=balance[winners]+winner_amount
bal[i, ]=balance
}
return(list(balance=bal,last_balance=balance))
}
set.seed(123)
hb_experiment(param)$last_balance
计算结果如下：
[1] 31.244667 82.693345 18.072119 44.561940 62.869733 33.399248 47.003107 45.550620
[9] 77.108569 70.435722 54.282843 26.981727 54.743774 80.302944 28.315862 43.980291
[17] 48.801944 82.694191 82.937560 -10.998868 34.301398 80.640606 60.678044 47.335579
[25] 40.126604 52.545811 23.386063 62.666160 92.201977 72.426845 41.550410 40.116285
[33] 50.507903 81.295672 51.165890 43.360551 34.928390 64.378636 42.700558 -8.903608
[41] 9.101540 78.605832 46.349353 64.175663 61.900195 13.606066 50.006777 68.509099
[49] 41.213098 54.141264
最后可以看出破产的玩家都是因为中头奖中的太多，导致入不敷出；相反，大赢家则是闷声发财这种。经统计，他获得第二名3次，第三名2次，第四名2次，第五名4次，等等。