英文:
Simulate an AR(1) process that approximates a gamma distribution
问题
我想模拟一个AR(1)时间序列，该序列近似于伽玛分布而不是正态分布。我希望结果具有指定的速率和形状，并具有给定的phi的AR(1)过程。我已经看到这个帖子，并理解严格的伽玛模拟存在问题，但我对近似形状和速率的结果会感到满意。粗略调整形状和速率，作为phi的函数，可能已经足够满足我的需求。有什么想法？
例如，
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'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.,
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)
[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 -> 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)
bar的ACF将与alpha更接近，而不是从随机生成的AR(1)中预期的：
# 随机高斯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
)
#> normal gamma
#> 1 1.000000000 1.00000000
#> 2 0.334034734 0.29929104
#> 3 0.060370711 0.09063563
#> 4 0.001948240 0.03818762
#> 5 -0.004354808 0.01196619
如果这种行为不理想，一个选项是将alpha设置为从随机生成的AR(1)的ACF（直到它不再单调下降）：
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
)
#> normal gamma gamma2
#> 1 1.000000000 1.00000000 1.000000000
#> 2 0.334034734 0.29929104 0.333494498
#> 3 0.060370711 0.09063563 0.060757048
#> 4 0.001948240 0.03818762 -0.017663549
#> 5 -0.004354808 0.01196619 0.008826011
为了方便，包括acf.reorder函数：
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))]
}
英文:
Borrowing the function from this CV answer, we can first sample from the desired distribution, then reorder the samples to approximate the expected ACF from the AR process.
set.seed(721893540)
n <- 500
foo <- rgamma(n = n, shape = 1.5 , rate = 5)
phi <- 0.3
# match the autocorrelation for the first 4 lags, when the
# autocorrelation (in the limit as n -> Inf) first drops below 0.01
alpha <- phi^(0:ceiling(log(0.01)/log(0.3)))
# reorder foo to get the desired ACF (ARIMA(0,0,1), phi = 0.3)
bar <- acf.reorder(foo, alpha)
The ACF of bar will match alpha more closely than would be expected from a randomly generated AR(1):
# autocorrelation function of a random Gaussian 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
)
#> normal gamma
#> 1 1.000000000 1.00000000
#> 2 0.334034734 0.29929104
#> 3 0.060370711 0.09063563
#> 4 0.001948240 0.03818762
#> 5 -0.004354808 0.01196619
If this behavior is undesirable, one option is to set alpha to the ACF from a randomly generated AR(1) (until it stops decreasing monotonically).
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
)
#> normal gamma gamma2
#> 1 1.000000000 1.00000000 1.000000000
#> 2 0.334034734 0.29929104 0.333494498
#> 3 0.060370711 0.09063563 0.060757048
#> 4 0.001948240 0.03818762 -0.017663549
#> 5 -0.004354808 0.01196619 0.008826011
Including acf.reorder for convenience.
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]*(w1[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))]
}
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论