Rstan - 高斯随机游走答案

【问题标题】：Rstan - Gaussian Random WalkRstan - 高斯随机游走
【发布时间】：2021-07-29 16:44:18
【问题描述】：

我正在尝试根据高斯随机游走对参数进行采样。在 R 中，代码如下所示：

#simulate a Gaussian random walk
#N  :  number of steps
#x0 :  initial offset
#mu :  drift velocity
#variance : step size
Gauss_RandomWalk <- function(N, x0, mu, variance) {
  z <- cumsum(rnorm(n=N, mean=mu, sd=sqrt(variance)))
  t <- 1:N
  x <- (x0 + t*mu + z)
  return(x)
}

实际上，设置x0=0.、mu=0.、variance=0.035**2时的结果看起来不错且合理：

[1]  0.040703269  0.009159686  0.052360030  0.059352074  0.092739218  0.098240752  0.113957813  0.064187776
[9]  0.062757728  0.063948224  0.034591074  0.004828493  0.019809969  0.032135111  0.025692763 -0.031678858
[17] -0.048033007 -0.020708105 -0.032231674  0.004917305  0.030961430  0.099054042  0.043441737 -0.010513085

每当我尝试在 Stan 中执行此操作时，例如根据 here 所代表的内容：

// model to be fitted
model {
  sqrtQ ~ student_t(2, 0, 0.035);
  // here we define the random walk for the log_Rt parameter
  log_Rt[1] ~ normal(0. , 0.035);
    for (t in 2:number_days) {
        log_Rt[t] ~ normal(log_Rt[t-1], sqrtQ);
  }
  print(log_Rt);
(...)
}

结果不是那么好。例如，跳跃-0.839043 -> 1.91956 大约是标准差的 85 倍，在统计上是不可能的……但是为什么会发生这种跳跃呢？

Chain 1: [0.382303,-0.489057,0.33374,-0.839043,1.91956,0.249953,-1.88793,1.61106,1.11189,-0.725063,-0.513174,1.79012,1.57758,0.75819,0.525524,1.44762,1.19118,0.485563,-1.48318,-1.36389,1.96355,0.321416,-0.365132,-0.644287,0.0981577,1.02943,-1.27993,-1.98085,-1.75191,-1.76489,1.60888,1.48925,-0.0452427,-0.92583,-1.21594,-0.906329,-0.700237,-0.208039,0.493656,0.490295,-1.61091,1.94587,0.758567,-1.02318,-1.92659,-0.999492,1.78042,0.214125,-0.0158054,-0.753422, .......

编辑：我也试过了：

// model to be fitted
model {
  // here we define the random walk for the log_Rt parameter
  log_Rt[1] ~ normal(0. , 0.035);
    for (t in 2:number_days) {
        log_Rt[t] ~ normal(log_Rt[t-1], 0.035);
  }
  print(log_Rt);
(...)
}

但事情根本没有改变。

【问题讨论】：

标签： r stan rstan

【解决方案1】：

您使用的正常函数实际上试图估计log_Rt 的概率，因此会产生奇怪的结果。您可以使用生成数量来执行类似于您在 r 中尝试的操作：

require(rstan)
#> Loading required package: rstan
#> Loading required package: StanHeaders
#> Loading required package: ggplot2
#> rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)
scode<-"
generated quantities {
 real y[10];
  y[1] = normal_rng(0, .035);
    for (t in 2:10) {
        y[t] = normal_rng(y[t-1], 0.035);
  }}"
m <- stan_model(model_code = scode) 
f<-sampling(m, iter=1, chain=1, algorithm='Fixed_param')
#> 
#> SAMPLING FOR MODEL '82a947a4909cd091541c3476a3eab1f0' NOW (CHAIN 1).
#> Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0 seconds (Warm-up)
#> Chain 1:                1.4e-05 seconds (Sampling)
#> Chain 1:                1.4e-05 seconds (Total)
#> Chain 1:
extract(f,'y')
#> $y
#>           
#> iterations       [,1]       [,2]       [,3]       [,4]         [,5]        [,6]
#>       [1,] 0.06017424 0.07683703 0.08551969 0.04428875 -0.001519853 0.000953893
#>           
#> iterations        [,7]        [,8]        [,9]       [,10]
#>       [1,] 0.008921459 -0.06060675 -0.04087786 -0.02217786

^{由reprex package (v2.0.1) 于 2021-09-28 创建}

【讨论】：