n <- 500
foo <- rgamma(n = n, shape = 1.5 , rate = 5)
phi <- 0.3
bar <- arima.sim(list(order = c(1,0,0), ar = phi), n = n, 
                 rand.gen = rgamma, shape = 1.5,
                 rate = 5)
ggplot() + 
  geom_density(aes(x=foo),fill="#FF8C00",alpha=0.5) + 
  geom_density(aes(x=bar),fill="#A034F0",alpha=0.5)

<details>
<summary>英文:</summary>
I’d like to simulate an AR(1) time series that approximates a gamma distribution rather than a normal distribution. I’d like the result to have a specified rate and shape along with a AR(1) process with a given phi. I have seen this [post][1] and understand that a strict gamma simulation is problematic, but I&#39;d be satisfied with a result that approximates the shape and rate. Some crude adjustment of shape and rate as a function of phi is likely enough for my needs. Any idea?
E.g.,

  [1]: https://stats.stackexchange.com/questions/180109/generate-a-random-variable-which-follow-gamma-distribution-and-ar1-process-sim
</details>
# 答案1
**得分**: 2
从[这个CV回答][1]中借用函数，我们可以首先从所需的分布中进行采样，然后 *重新排序* 这些样本以近似来自AR过程的期望ACF。
```R
set.seed(721893540)
n <- 500
foo <- rgamma(n = n, shape = 1.5 , rate = 5)
phi <- 0.3
# 匹配前4个滞后的自相关，当自相关（当n -&gt; Inf时）首次下降到0.01以下
alpha <- phi^(0:ceiling(log(0.01)/log(0.3)))
# 重新排序foo以获得所需的ACF（ARIMA(0,0,1)，phi = 0.3）
bar <- acf.reorder(foo, alpha)

# 随机高斯AR(1)的自相关函数：
ACF <- acf(arima.sim(list(order = c(1,0,0), ar = phi), n), plot = 0)$acf
data.frame(
  normal = ACF[1:length(alpha)],
  gamma = acf(bar, length(alpha) - 1, plot = 0)$acf
)
#&gt;         normal      gamma
#&gt; 1  1.000000000 1.00000000
#&gt; 2  0.334034734 0.29929104
#&gt; 3  0.060370711 0.09063563
#&gt; 4  0.001948240 0.03818762
#&gt; 5 -0.004354808 0.01196619

alpha <- ACF[1:which.max(diff(ACF) > 0)]
bar2 <- acf.reorder(foo, alpha)
data.frame(
  normal = ACF[1:length(alpha)],
  gamma = acf(bar, length(alpha) - 1, plot = 0)$acf,
  gamma2 = acf(bar2, length(alpha) - 1, plot = 0)$acf
)
#&gt;         normal      gamma       gamma2
#&gt; 1  1.000000000 1.00000000  1.000000000
#&gt; 2  0.334034734 0.29929104  0.333494498
#&gt; 3  0.060370711 0.09063563  0.060757048
#&gt; 4  0.001948240 0.03818762 -0.017663549
#&gt; 5 -0.004354808 0.01196619  0.008826011

acf.reorder <- function(x, alpha) {
  tol <- 1e-5
  maxIter <- 10L
  n <- length(x)
  xx <- sort(x)
  y <- rnorm(n)
  w0 <- w <- alpha1 <- alpha
  m <- length(alpha)
  i1 <- sequence((m - 1):1)
  i2 <- sequence((m - 1):1, 2:m)
  i3 <- cumsum((m - 1):1)
  tol10 <- tol/10
  iter <- 0L
  x <- xx[rank(filter(y, w, circular = TRUE))]
  SSE0 <- Inf
  f <- function(ww) {
    sum((c(1, diff(c(0, cumsum(ww[i1]*(ww[i2]))[i3]))/sum(ww^2)) - alpha1)^2)
  }
  ACF <- function(x) acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m]
  
  while ((SSE <- sum((ACF(x) - alpha)^2)) > tol) {
    if (SSE < SSE0) {
      SSE0 <- SSE
      w <- w0
    }
    if ((iter <- iter + 1L) == maxIter) break
    w1 <- w0
    a <- 0
    sse0 <- Inf
    
    while (max(abs(alpha1 - a)) > tol10) {
      a <- c(1, diff(c(0, cumsum(w1[i1]*(ww1[i2]))[i3]))/sum(w1^2))
      
      if ((sse <- sum((a - alpha1)^2)) < sse0) {
        sse0 <- sse
        w0 <- w1
      } else {
        # w0 failed to converge; try optim
        w1 <- optim(w0, f, method = "L-BFGS-B")$par
        a <- c(1, diff(c(0, cumsum(w1[i1]*(w1[i2]))[i3]))/sum(w1^2))
        if (sum((a - alpha1)^2) < sse0) w0 <- w1
        break
      }
      
      w1 <- (w1*alpha1/a + w1)/2
    }
    
    x <- xx[rank(filter(y, w0, circular = TRUE))]
    alpha1 <- (alpha1*alpha/ACF(x) + alpha1)/2
  }
  
  xx[rank(filter(y, w, circular = TRUE))]
}