在R中模拟中式餐厅流程

Question

我正在尝试在R中模拟Chinese Restaurant process，并且想知道我是否可以对这个粗糙的实现进行任何效率改进。

iTables = 200  # number of tables
iSampleSize = 1000  # number of diners

# initialize the list of tables
listTableOccupants = vector('list', iTables)

for(currentDiner in seq.int(iSampleSize)) {
  # occupation probabilities for the next diner
  vProbabilities = sapply(listTableOccupants, 
                          function(x) ifelse(!is.null(x),
                                             length(x)/currentDiner,
                                             1/currentDiner))
  # pick the index of the lucky table
  iTable = sample.int(iTables, size = 1, prob = vProbabilities)

  # add to the list element corresponding to the table
  listTableOccupants[[iTable]] = 
    c(listTableOccupants[[iTable]], currentDiner) 
}

我特别关注的是这一行：

  # add to the list element corresponding to the table
  listTableOccupants[[iTable]] = 
    c(listTableOccupants[[iTable]], currentDiner) 

这是有效的吗？

Answer 1

为了避免空间重新分配和稀疏数据结构，您可以改为对每个用餐者应用一个表标签。例如,

nDnr <- 100  # number of diners; must be at least 2
vDnrTbl <- rep(0, nDnr)  # table label for each diner
alpha <- 2  # CRP parameter

vDnrTbl[1] <- 1
for (dnr in 2:length(vDnrTbl)) {

    # compute occupation probabilities for current diner
    vOcc <- table(vDnrTbl[1:(dnr-1)])
    vProb <- c(vOcc, alpha) / (dnr - 1 + alpha)

    # add table label to diner
    nTbl <- as.numeric(names(vOcc)[length(vOcc)])  # avoid overhead of finding max of possibly large vector
    vDnrTbl[dnr] <- sample.int(nTbl+1, size=1, prob=vProb)
}

您可以从vDnrTbl获取listTableOccupants

nTbl <- max(c(nTbl, vDnrTbl[dnr]))
listTableOccupants <- lapply(1:nTbl, function(t) which(vDnrTbl == t))