模拟数据采样的笔记

在计算机模拟中,往往需要对一条轨迹中得到的数据求平均值.为了证明结论可靠,给出这个平均值的误差估计是十分必要的.由于模拟数据具有很高的自相关性,误差估计是一个并不简便的计算.本文首先通过一个简单的玩具模型构建一套自相关的数据,以此为基础介绍了平均值的误差估计的理论和计算,并讨论了如何在实际工作中选择最佳的采样数目.

1.数据生成

Markov Chain Monte Carlo(MCMC) 算法

考虑一个两态系统, 分别用0/1表示,其转移矩阵可以用一个2x2的跃迁矩阵表示 \[ \begin{bmatrix} 1-p & p \\ p & 1-p \end{bmatrix} \] 其中p为系统在下一步中离开原来状态的概率.

系统依据跃迁矩阵变化的代码如下:

simula <- function(trans, N) {
    transita <- function(char, trans) {
        sample(colnames(trans), 1, prob = trans[char, ])
    }
    sim <- integer(N)
    sim[1] <- as.integer(sample(colnames(trans), 1))
    for (i in 2:N) {
        sim[i] <- as.integer(transita(as.character(sim[i - 1]), trans))
    }
    return(sim)
}

生成数据集

下文中,我们以p=0.1为例,生成一个包含10000个采样的数据集.

tmax <- 10000
p <- 0.1
mat <- matrix(c(1 - p, p, p, 1 - p), ncol = 2, byrow = TRUE)
colnames(mat) <- c("0", "1")
row.names(mat) <- c("0", "1")
data <- simula(mat, tmax)

生成的数据为0和1组成的数组,如

print(data[1:70])

##  [1] 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [36] 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

2.本模型参数的解析解

通过简单计算我们可以知道,数据集的期望为 \[ \left \langle x \right \rangle = (\frac {N}{2} \times 0 + \frac{N}{2} \times 1) / N = 0.5 \] 方差为 \[ \left \langle x^2 \right \rangle - \left \langle x \right \rangle ^2 = 0.5 - (0.5)^2 = 0.25 \] 在计算机模拟等领域,常见的一个错误是将方差作为期望的误差估计给出了.实际上,很多时候当我们得到一个物理量,也就是对一个均值进行了估计之后,关心的是这个估计的精度,也就是这个平均值的误差,而非原随机变量的误差.对于一个涨落很大的随机变量(如压强),其平均值或许已经估计得很准,远远小于其涨落大小,那究竟这个精度该如何计算见后文.

3.数据集的采样性质

首先给出采样冗余的计算公式(参考文献 Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations, by Chodera, JCTC, 3-26, 2007)

calstateff <- function(skip = 1, N = 10, plt = FALSE) {
    ac <- acf(data[(1:(tmax/skip)) * skip], lag.max = lmax/skip, plot = plt)
    tau <- sum((1 - (ac$lag/lmax)) * ac$acf) - 1
    if (tau < 0) 
        tau <- 0  #由于数据的噪音有时tau可能取到负值
    gf <- 1 + 2 * tau
    print(sprintf("skip=%d, Neff=%f", skip, gf))
    return(gf)
}

我们取最大时间范围为lmax=400,当取全部数据进行计算时

lmax <- 400
gfdat <- numeric(20)
gfdat[1] <- calstateff(1, lmax, TRUE)

plot of chunk b0

## [1] "skip=1, Neff=10.721383"

可以看到数据冗余量约为10左右,采样效率很低.逐渐加大采样间隔可以得到

for (skip in 2:20) {
    gfdat[skip] <- calstateff(skip, lmax)
}

## [1] "skip=2, Neff=5.681595"
## [1] "skip=3, Neff=4.269735"
## [1] "skip=4, Neff=3.360553"
## [1] "skip=5, Neff=2.432393"
## [1] "skip=6, Neff=2.277844"
## [1] "skip=7, Neff=1.816193"
## [1] "skip=8, Neff=1.669858"
## [1] "skip=9, Neff=2.078204"
## [1] "skip=10, Neff=1.598440"
## [1] "skip=11, Neff=1.568251"
## [1] "skip=12, Neff=1.714554"
## [1] "skip=13, Neff=1.532862"
## [1] "skip=14, Neff=1.000000"
## [1] "skip=15, Neff=1.354868"
## [1] "skip=16, Neff=1.306202"
## [1] "skip=17, Neff=1.003490"
## [1] "skip=18, Neff=1.590329"
## [1] "skip=19, Neff=1.625481"
## [1] "skip=20, Neff=1.180311"

在skip大于10后,采样效率达到最大,冗余接近1即没有冗余,采样效率达到最高.

4.Bootstrapping方法

上一节中我们看到,对于一组数据由于其自身的相关性,并不需要取所有的数据点进行计算即可保证信息不会丢失(也就是精度不变),更多的采样并不能增加精度,而过少的采样会造成精度的下降.这一节,我们用bootstrapping方法来验证这一结果.

有放回的重抽样

对于初始的数据集合,可以通过bootstrap方法来估计其平均值的精度.本方法的核心思想在于,以相同概率从原数据集中抽取采样,并将其再放回去继续这一过程,得到一组新的数据,多次重复统计所得均值的分布,即可描述原数据集的性质.代码如下

bootstrapping <- function(dat, n = 100) {
    m <- numeric(n)
    for (i in 1:n) {
        m[i] <- mean(sample(dat, replace = TRUE))
    }
    return(m)
}

我们以不同的采样间隔构建数据集合,并评估其均值的精度

vardat <- numeric(20)
for (skip in 1:20) {
    subdat <- data[(1:(tmax/skip)) * skip]
    re <- bootstrapping(subdat, 500)
    vardat[skip] <- var(re)
    print(sprintf("mean_raw=%f,mean_re=%f, var_re=%f", mean(data), mean(re), 
        vardat[skip]))
}

## [1] "mean_raw=0.475400,mean_re=0.475098, var_re=0.000024"
## [1] "mean_raw=0.475400,mean_re=0.473792, var_re=0.000054"
## [1] "mean_raw=0.475400,mean_re=0.474940, var_re=0.000073"
## [1] "mean_raw=0.475400,mean_re=0.473638, var_re=0.000105"
## [1] "mean_raw=0.475400,mean_re=0.470961, var_re=0.000118"
## [1] "mean_raw=0.475400,mean_re=0.476696, var_re=0.000150"
## [1] "mean_raw=0.475400,mean_re=0.481789, var_re=0.000175"
## [1] "mean_raw=0.475400,mean_re=0.469251, var_re=0.000193"
## [1] "mean_raw=0.475400,mean_re=0.470734, var_re=0.000224"
## [1] "mean_raw=0.475400,mean_re=0.461850, var_re=0.000250"
## [1] "mean_raw=0.475400,mean_re=0.484262, var_re=0.000256"
## [1] "mean_raw=0.475400,mean_re=0.476992, var_re=0.000320"
## [1] "mean_raw=0.475400,mean_re=0.479347, var_re=0.000297"
## [1] "mean_raw=0.475400,mean_re=0.468571, var_re=0.000360"
## [1] "mean_raw=0.475400,mean_re=0.487850, var_re=0.000386"
## [1] "mean_raw=0.475400,mean_re=0.484570, var_re=0.000394"
## [1] "mean_raw=0.475400,mean_re=0.471850, var_re=0.000428"
## [1] "mean_raw=0.475400,mean_re=0.466984, var_re=0.000425"
## [1] "mean_raw=0.475400,mean_re=0.470388, var_re=0.000489"
## [1] "mean_raw=0.475400,mean_re=0.455916, var_re=0.000528"

实际的精度,还需要乘以采样冗余,我们给出精度随采样间隔的变化

plot(1:20, gfdat * vardat)

plot of chunk unnamed-chunk-7

虽然有较大噪音,但我们依然可以看出,在skip<10左右的区域,由于没有信息丢失,均值的误差处在同一水平.而skip继续增大时,误差线性增加,也就发生了信息的丢失.