感谢@jblood94提供了他对上一个问题的答案的链接，这个问题使用这里提供的acf.reorder函数提供了一个非常好的解决方案。
我已经调整并简化了下面的函数，以生成从给定的负二项分布（由mu和size指定）绘制的N自相关计数的时间序列。

acf.negbin <- function(N, mu, size, alpha, max.iter = 100, tol = 1e-5) {
  
  m = length(alpha)
  
  generate = function(){
   x = sort(rnbinom(N,size=size,mu=mu))
   y <- rnorm(length(x))
   x[rank(stats::filter(y, alpha, circular = TRUE))]
  }
  
  a = generate()
  iter <- 0L
  
  ACF <- function(x) acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m]
  SSE <- function(x,alpha) sum((ACF(x) - alpha)^2)
    
  while ((SSE0 <- SSE(a, alpha)) > tol) {
    if ((iter <- iter + 1L) == max.iter) break
    a1 <- generate()
    if(SSE(a1,alpha) < SSE0) a <- a1
  }
    
  return(a)
  
}

如原始帖子中所述，该方法的工作原理是重复从负二项分布中抽取样本，并对它们进行重新排序，以近似于alpha给出的所需自相关函数。在修改后的函数中，我已经将随机计数生成移动到while循环中，这似乎允许它更有效地为短时间序列找到一个很好的近似。可以通过增加迭代次数（max.iter）或减小容差（tol）来获得更好的近似值，容差是期望值alpha和观测值alpha之间的平方误差之和，在此值下可以接受解。
测试表明，它非常接近我上面的示例所需的时间序列

set.seed(2023)
x = acf.negbin(25,mu=10000,size=5.217,alpha=c(1,-0.5),max.iter=1000)
plot(x,type='l')
mean(x)
acf(x)$acf[1:2]

再次感谢@jblood94的解决方案，并为每个原始帖子+1！